##########################################################################################
# Script for the "Sources of Public Opinion (Mis)perceptions by Elected Representatives" #
# By: Simon Otjes                                                                        #
# Date 30/8/2024                                                                         #
# Version 5.0                                                                            #
##########################################################################################

# Hi! This code comes in 5 parts
# PART 1 reads in packages & data
# PART 2 makes the variables
# PART 3 makes the descriptives
# PART 4 performs the regression analyses
# PART 5 performs the mediation analyses

########################################################################################
# PART 1: READING IN
# PART 1.1: PACKAGES & FUNCTIONS

rm(list=ls())

library(mokken)
library(lme4)
library(texreg)
library(mediation)
library(performance)
library(glmmTMB)

# a function generating summary data
overview<-function(vector) {

	return<-c(round(mean(vector,na.rm=TRUE),2),
	round(median(vector,na.rm=TRUE),2),
	round(sd(vector,na.rm=TRUE),2),
	round(min(vector,na.rm=TRUE),2),
	round(max(vector,na.rm=TRUE),2),
	length(vector)-sum(is.na(vector)))
	return(return)
}

# set the seed (relevant for the mediation)
set.seed(100)

########################################################################################
# PART 1.2: ELITE-LEVEL DATA

# set working directory
setwd("~/Dropbox/Data Anne Simon/Perceptions paper/Replication file/")

# read in responses
DATA01<-read.csv("datasets/240529 dataset pseudonymous.csv")

########################################################################################
# Part 1.3 RESPONSE RATE

# Sample is the document of all the people we sent an email to

fullsample<-read.csv("datasets/240529 full sample pseudonyms.csv")

#remove double entries

sample<-fullsample[fullsample$Doublures==0,]

# obtain unique IDs of respondents

runique<-(DATA01$uniqueid)
runi<-cbind(as.vector(runique),1:length(runique),1)
colnames(runi)<-c("UniqueID","No","Participated")

# merge two files

rsample<-merge(x=sample,y=runi,all.x=TRUE,all.y=TRUE,by="UniqueID")

# respondents who have missing data for participation have not participated
rsample$Participated[is.na(rsample$Participated)]<-0

# generate Table A1 on response rates
# emails sent (Table A1/col 1)
xtabs(~rsample$Sample)

sum(xtabs(~rsample$Sample)[1:3])
sum(xtabs(~rsample$Sample)[4:6])

sum(xtabs(~rsample$Sample)[c(1,4)])
sum(xtabs(~rsample$Sample)[c(2,5)])
sum(xtabs(~rsample$Sample)[c(3,6)])

sum(xtabs(~rsample$Sample))

#identify bad emails (e.g. bounces)
actual<-rsample$In.Survey
actual[as.logical(rsample$Email.missing)]<-0

#actual (Table A1/col 2)
xtabs(~rsample$Sample[actual==1])

sum(xtabs(~rsample$Sample[actual==1])[1:3])
sum(xtabs(~rsample$Sample[actual==1])[4:6])

sum(xtabs(~rsample$Sample[actual==1])[c(1,4)])
sum(xtabs(~rsample$Sample[actual==1])[c(2,5)])
sum(xtabs(~rsample$Sample[actual==1])[c(3,6)])

sum(xtabs(~rsample$Sample[actual==1]))

# responses (Table A1/col 3)

xtabs(~rsample$Sample[rsample$Participated==1])

sum(xtabs(~rsample$Sample[rsample$Participated==1])[1:3])
sum(xtabs(~rsample$Sample[rsample$Participated==1])[4:6])

sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(1,4)])
sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(2,5)])
sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(3,6)])

sum(xtabs(~rsample$Sample[rsample$Participated==1]))

# updated sample (Table A1/col 4)

updated<-as.numeric((rsample$Participated==1)|(actual==1))

xtabs(~rsample$Sample[updated==1])

sum(xtabs(~rsample$Sample[updated==1])[1:3])
sum(xtabs(~rsample$Sample[updated==1])[4:6])

sum(xtabs(~rsample$Sample[updated==1])[c(1,4)])
sum(xtabs(~rsample$Sample[updated==1])[c(2,5)])
sum(xtabs(~rsample$Sample[updated==1])[c(3,6)])

sum(xtabs(~rsample$Sample[updated==1]))

#response rate (Table A1/col 5)
round(xtabs(~rsample$Sample[rsample$Participated==1])/xtabs(~rsample$Sample[updated==1])*100,1)

round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[1:3])/sum(xtabs(~rsample$Sample[updated==1])[1:3])*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[4:6])/sum(xtabs(~rsample$Sample[updated==1])[4:6])*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(1,4)])/sum(xtabs(~rsample$Sample[updated==1])[c(1,4)])*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(2,5)])/sum(xtabs(~rsample$Sample[updated==1])[c(2,5)])*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(3,6)])/sum(xtabs(~rsample$Sample[updated==1])[c(3,6)])*100,1)

round(sum(xtabs(~rsample$Sample[rsample$Participated==1]))/sum(xtabs(~rsample$Sample[updated==1]))*100,1)

#share of sample (Table A1/col 6)
round(xtabs(~rsample$Sample[rsample$Participated==1])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)

round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[1:3])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[4:6])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)

round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(1,4)])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(2,5)])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(3,6)])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)
round(sum(xtabs(~rsample$Sample[rsample$Participated==1])[c(1:6)])/sum(xtabs(~rsample$Sample[rsample$Participated==1]))*100,1)

# generate Table A2 on representativity
# share of population (Table A2/col 1) for party
round(xtabs(~rsample$Party.1)/sum(xtabs(~rsample$Party.1))*100,0)

# to check values below 1: this should be rounded as "<1"
round(xtabs(~rsample$Party.1)/sum(xtabs(~rsample$Party.1))*100,1)[8]

#share of respondents (Table A2/col 2)
round(xtabs(~rsample$Party.1[rsample$Participated==1])/sum(xtabs(~rsample$Party.1[rsample$Participated==1]))*100,0)

# to check values below 1: this should be rounded as "<1"
round(xtabs(~rsample$Party.1[rsample$Participated==1])/sum(xtabs(~rsample$Party.1[rsample$Participated==1]))*100,1)[8]

#share of population (Table A2/col 1) for gender
round(xtabs(~rsample$morf)/sum(xtabs(~rsample$morf))*100,0)

#share of respondents (Table A2/col 2) for gender
round(xtabs(~rsample$morf[rsample$Participated==1])/sum(xtabs(~rsample$morf[rsample$Participated==1]))*100,0)

########################################################################################
# generate data set on the basis of survey answers and contextual variables in file with all respondents
colnames(DATA01)[1]<-"UniqueID"

DATA<-merge(x=DATA01,y=sample,by="UniqueID",all.x=TRUE,all.y=FALSE)

########################################################################################
# PART 1.4: VOTER-LEVEL DATA

DK<-read.csv("datasets/240529 mass sample denmark.csv")
NL<-read.csv("datasets/240529 mass sample netherlands.csv")

########################################################################################
# PART 2: MAKING VARIABLES
# PART 2.2: DEPENDENT VARIABLES
# there are six DVs: bias, accuracy and a dichotomous measure for weighted and unweighted data

# estimates of politicians

EstASYL<-DATA$NationalPosition.Asylumseekers
EstCCTV<-DATA$NationalPosition.CCTV
EstCULT<-DATA$NationalPosition.Culture
EstCARE<-DATA$NationalPosition.ElderlyCare
EstROAD<-DATA$NationalPosition.TrafficJams

# actual shares

TrwASYL<-rep(NA,dim(DATA)[1])
TrwCCTV<-rep(NA,dim(DATA)[1])
TrwCULT<-rep(NA,dim(DATA)[1])
TrwCARE<-rep(NA,dim(DATA)[1])
TrwROAD<-rep(NA,dim(DATA)[1])

# Question order in mass surveys
# asyl = 6
# cctv = 10
# cult = 8
# care = 2
# road = 1

# calculate the shares answering positively ("Grotendeels/Helemaal mee eens" in Dutch and ("Helt/Stort set enig" in Danish) while removing don't knows ("Weet ik niet" in Dutch and "Ved ikke" in Danish)
TrwASYL[DATA$Country=="NL"]<-(sum(NL$v7a_6=="Grotendeels mee eens")+sum(NL$v7a_6=="Helemaal mee eens"))/(length(NL$v7a_6)-sum(NL$v7a_6=="Weet ik niet"))*100
TrwASYL[DATA$Country=="DK"]<-(sum(DK$q1_6=="Helt enig")+sum(DK$q1_6=="Stort set enig"))/(length(DK$q1_6)-sum(DK$q1_6=="Ved ikke"))*100

TrwCCTV[DATA$Country=="NL"]<-(sum(NL$v7a_10=="Grotendeels mee eens")+sum(NL$v7a_10=="Helemaal mee eens"))/(length(NL$v7a_10)-sum(NL$v7a_10=="Weet ik niet"))*100
TrwCCTV[DATA$Country=="DK"]<-(sum(DK$q1_10=="Helt enig")+sum(DK$q1_10=="Stort set enig"))/(length(DK$q1_10)-sum(DK$q1_10=="Ved ikke"))*100

TrwCULT[DATA$Country=="NL"]<-(sum(NL$v7a_8=="Grotendeels mee eens")+sum(NL$v7a_8=="Helemaal mee eens"))/(length(NL$v7a_8)-sum(NL$v7a_8=="Weet ik niet"))*100
TrwCULT[DATA$Country=="DK"]<-(sum(DK$q1_8=="Helt enig")+sum(DK$q1_8=="Stort set enig"))/(length(DK$q1_8)-sum(DK$q1_8=="Ved ikke"))*100

TrwCARE[DATA$Country=="NL"]<-(sum(NL$v7a_2=="Grotendeels mee eens")+sum(NL$v7a_2=="Helemaal mee eens"))/(length(NL$v7a_2)-sum(NL$v7a_2=="Weet ik niet"))*100
TrwCARE[DATA$Country=="DK"]<-(sum(DK$q1_2=="Helt enig")+sum(DK$q1_2=="Stort set enig"))/(length(DK$q1_2)-sum(DK$q1_2=="Ved ikke"))*100

TrwROAD[DATA$Country=="NL"]<-(sum(NL$v7a_1=="Grotendeels mee eens")+sum(NL$v7a_1=="Helemaal mee eens"))/(length(NL$v7a_1)-sum(NL$v7a_1=="Weet ik niet"))*100
TrwROAD[DATA$Country=="DK"]<-(sum(DK$q1_1=="Helt enig")+sum(DK$q1_1=="Stort set enig"))/(length(DK$q1_1)-sum(DK$q1_1=="Ved ikke"))*100

# actual shares for weighted data

TruASYL<-rep(NA,dim(DATA)[1])
TruCCTV<-rep(NA,dim(DATA)[1])
TruCULT<-rep(NA,dim(DATA)[1])
TruCARE<-rep(NA,dim(DATA)[1])
TruROAD<-rep(NA,dim(DATA)[1])

# calculate the shares answering positively while removing don't knows and applying weights

TruASYL[DATA$Country=="NL"]<-(sum((NL$v7a_6=="Grotendeels mee eens")*NL$weging_totaal_pol)+sum((NL$v7a_6=="Helemaal mee eens")*NL$weging_totaal_pol))/(sum(NL$weging_totaal_pol)-sum((NL$v7a_6=="Weet ik niet")*NL$weging_totaal_pol))*100
TruASYL[DATA$Country=="DK"]<-(sum((DK$q1_6=="Helt enig")*DK$weight)+sum((DK$q1_6=="Stort set enig")*DK$weight))/(sum(DK$weight)-sum((DK$q1_6=="Ved ikke")*DK$weight))*100

TruCCTV[DATA$Country=="NL"]<-(sum((NL$v7a_10=="Grotendeels mee eens")*NL$weging_totaal_pol)+sum((NL$v7a_10=="Helemaal mee eens")*NL$weging_totaal_pol))/(sum(NL$weging_totaal_pol)-sum((NL$v7a_10=="Weet ik niet")*NL$weging_totaal_pol))*100
TruCCTV[DATA$Country=="DK"]<-(sum((DK$q1_10=="Helt enig")*DK$weight)+sum((DK$q1_10=="Stort set enig")*DK$weight))/(sum(DK$weight)-sum((DK$q1_10=="Ved ikke")*DK$weight))*100

TruCULT[DATA$Country=="NL"]<-(sum((NL$v7a_8=="Grotendeels mee eens")*NL$weging_totaal_pol)+sum((NL$v7a_8=="Helemaal mee eens")*NL$weging_totaal_pol))/(sum(NL$weging_totaal_pol)-sum((NL$v7a_8=="Weet ik niet")*NL$weging_totaal_pol))*100
TruCULT[DATA$Country=="DK"]<-(sum((DK$q1_8=="Helt enig")*DK$weight)+sum((DK$q1_8=="Stort set enig")*DK$weight))/(sum(DK$weight)-sum((DK$q1_8=="Ved ikke")*DK$weight))*100

TruCARE[DATA$Country=="NL"]<-(sum((NL$v7a_2=="Grotendeels mee eens")*NL$weging_totaal_pol)+sum((NL$v7a_2=="Helemaal mee eens")*NL$weging_totaal_pol))/(sum(NL$weging_totaal_pol)-sum((NL$v7a_2=="Weet ik niet")*NL$weging_totaal_pol))*100
TruCARE[DATA$Country=="DK"]<-(sum((DK$q1_2=="Helt enig")*DK$weight)+sum((DK$q1_2=="Stort set enig")*DK$weight))/(sum(DK$weight)-sum((DK$q1_2=="Ved ikke")*DK$weight))*100

TruROAD[DATA$Country=="NL"]<-(sum((NL$v7a_1=="Grotendeels mee eens")*NL$weging_totaal_pol)+sum((NL$v7a_1=="Helemaal mee eens")*NL$weging_totaal_pol))/(sum(NL$weging_totaal_pol)-sum((NL$v7a_1=="Weet ik niet")*NL$weging_totaal_pol))*100
TruROAD[DATA$Country=="DK"]<-(sum((DK$q1_1=="Helt enig")*DK$weight)+sum((DK$q1_1=="Stort set enig")*DK$weight))/(sum(DK$weight)-sum((DK$q1_1=="Ved ikke")*DK$weight))*100

# 2.2.1 dichtomous approach
# dichtomous approach for weighted data

DicASYL<-(round(TruASYL/50,0)-1) == (round(EstASYL/50,0)-1) 
DicASYL[EstASYL==50] <- TRUE
DicCCTV<-(round(TruCCTV/50,0)-1) == (round(EstCCTV/50,0)-1) 
DicCCTV[EstCCTV==50] <- TRUE
DicCULT<-(round(TruCULT/50,0)-1) == (round(EstCULT/50,0)-1) 
DicCULT[EstCULT==50] <- TRUE
DicCARE<-(round(TruCARE/50,0)-1) == (round(EstCARE/50,0)-1) 
DicCARE[EstCARE==50] <- TRUE
DicROAD<-(round(TruROAD/50,0)-1) == (round(EstROAD/50,0)-1) 
DicROAD[EstROAD==50] <- TRUE

Dic<-c(DicASYL,DicCCTV,DicCULT,DicCARE,DicROAD)

# dichtomous approach for unweighted data

DwcASYL<-(round(TrwASYL/50,0)-1) == (round(EstASYL/50,0)-1) 
DwcASYL[EstASYL==50] <- TRUE
DwcCCTV<-(round(TrwCCTV/50,0)-1) == (round(EstCCTV/50,0)-1) 
DwcCCTV[EstCCTV==50] <- TRUE
DwcCULT<-(round(TrwCULT/50,0)-1) == (round(EstCULT/50,0)-1) 
DwcCULT[EstCULT==50] <- TRUE
DwcCARE<-(round(TrwCARE/50,0)-1) == (round(EstCARE/50,0)-1) 
DwcCARE[EstCARE==50] <- TRUE
DwcROAD<-(round(TrwROAD/50,0)-1) == (round(EstROAD/50,0)-1) 
DwcROAD[EstROAD==50] <- TRUE

Dwc<-c(DwcASYL,DwcCCTV,DwcCULT,DwcCARE,DwcROAD)

# 2.2.2 Bias
# Difference for weighted data

DifASYL<-EstASYL-TruASYL
DifCCTV<-EstCCTV-TruCCTV
DifCULT<-EstCULT-TruCULT
DifCARE<-EstCARE-TruCARE
DifROAD<-EstROAD-TruROAD

# Difference for unweighted data

DwfASYL<-EstASYL-TrwASYL
DwfCCTV<-EstCCTV-TrwCCTV
DwfCULT<-EstCULT-TrwCULT
DwfCARE<-EstCARE-TrwCARE
DwfROAD<-EstROAD-TrwROAD

# Flipping vectors so that they are all in left-right direction for weighted data

RightROAD<-  DifROAD
RightCARE<- -DifCARE
RightASYL<- -DifASYL
RightCCTV<-  DifCCTV
RightCULT<-  DifCULT

RightDif<-c(RightASYL,RightCCTV,RightCULT,RightCARE,RightROAD)

# Flipped vectors so that they are all in left-right direction for unweighted data

RwghtROAD<-  DwfROAD
RwghtCARE<- -DwfCARE
RwghtASYL<- -DwfASYL
RwghtCCTV<-  DwfCCTV
RwghtCULT<-  DwfCULT

RwghtDif<-c(RwghtASYL,RwghtCCTV,RwghtCULT,RwghtCARE,RwghtROAD)

# 2.2.3 Accuracy
# Absolute difference for weighted data

AbsASYL<-100-abs(DifASYL)
AbsCCTV<-100-abs(DifCCTV)
AbsCULT<-100-abs(DifCULT)
AbsCARE<-100-abs(DifCARE)
AbsROAD<-100-abs(DifROAD)

# Absolute difference for unweighted data

AwsASYL<-100-abs(DwfASYL)
AwsCCTV<-100-abs(DwfCCTV)
AwsCULT<-100-abs(DwfCULT)
AwsCARE<-100-abs(DwfCARE)
AwsROAD<-100-abs(DwfROAD)

Abs<-c(AbsASYL,AbsCCTV,AbsCULT,AbsCARE,AbsROAD)
Aws<-c(AwsASYL,AwsCCTV,AwsCULT,AwsCARE,AwsROAD)

# transformed accuracy
Sba<-(100-Abs)^(1/2)
Swa<-(100-Aws)^(1/2)

# divided by 100 to aid interpretation in some analyses
Abs1<-Abs/100
Aws1<-Aws/100

########################################################################################
# PART 2.3: GROUP VARIABLES
# There are nine group based variables: business and issue-specific societal engagement, business/firm, issue-specific societal, broad and narrow societal contacts, business and NGO information and the difference between the last two

# 2.3.1   Engagement
# 2.3.1.1 Business group engagement
linksBusiness<-cbind(	DATA$Links.BusinessGroups.Member,
						DATA$Links.BusinessGroups.Activities,
						DATA$Links.BusinessGroups.Donated,
						DATA$Links.BusinessGroups.VolunteerWork,
						DATA$Links.BusinessGroups.PaidWork,
						DATA$Links.BusinessGroups.None,
						DATA$Links.BusinessGroups.Don.t.know)

linksBusiness2<-linksBusiness[,1:5]
linksBusiness3<-linksBusiness
linksBusiness3[is.na(linksBusiness3)]<-0

# sixth column = "none", then the answer should be none
linksBusiness2[linksBusiness3[,6]==1,]<-0

# seventh column = "don't know", then the answer should be don't know
linksBusiness2[linksBusiness3[,7]==1,]<-NA
linksBusiness2[is.na(linksBusiness[,1:5])]<-NA

BusinessLinks<-rowMeans(linksBusiness2)

BL<-c(BusinessLinks,BusinessLinks,BusinessLinks,BusinessLinks,BusinessLinks)

# Create a binary variable

BusinessBinary<-1-as.numeric(BusinessLinks==0)
BusinessBinary[is.na(BusinessLinks)]<-NA

BB<-c(BusinessBinary,BusinessBinary,BusinessBinary,BusinessBinary,BusinessBinary)

# 2.3.1.2 Issue-specific, societal engagement
# the structure is the same as above with assigning none to 0, don't know to NA.
# human rights group engagement

linksHURI<-cbind(		DATA$Links.HumanRights.Member,
						DATA$Links.HumanRights.Activities,
						DATA$Links.HumanRights.Donated,
						DATA$Links.HumanRights.VolunteerWork,
						DATA$Links.HumanRights.PaidWork,
						DATA$Links.HumanRights.None,
						DATA$Links.HumanRights.Don.t.know)

linksHURI5<-linksHURI[,1:5]

# sixth column = "none"
# seventh column = "don't know"
LH6<-linksHURI[,6]==1
LH6[is.na(LH6)]<-FALSE

LH7<-linksHURI[,7]==1
LH7[is.na(LH7)]<-FALSE

linksHURI5[LH6,]<-0
linksHURI5[LH7,]<-NA

# culture group engagement

linksCULT<-cbind(		DATA$Links.Culture.Member,
						DATA$Links.Culture.Activities,
						DATA$Links.Culture.Donated,
						DATA$Links.Culture.VolunteerWork,
						DATA$Links.Culture.PaidWork,
						DATA$Links.Culture.None,
						DATA$Links.Culture.Don.t.know)

linksCULT5<-linksCULT[,1:5]

# sixth column = "none"
# seventh column = "don't know"

LC6<-linksCULT[,6]==1
LC6[is.na(LC6)]<-FALSE

LC7<-linksCULT[,7]==1
LC7[is.na(LC7)]<-FALSE

linksCULT5[LC6,]<-0
linksCULT5[LC7,]<-NA

# seniors' group engagement

linksSENI<-cbind(		DATA$Links.Seniors.Member,
						DATA$Links.Seniors.Activities,
						DATA$Links.Seniors.Donated,
						DATA$Links.Seniors.VolunteerWork,
						DATA$Links.Seniors.PaidWork,
						DATA$Links.Seniors.None,
						DATA$Links.Seniors.Don.t.know)

linksSENI5<-linksSENI[,1:5]

# sixth column = "none"
# seventh column = "don't know"

LS6<-linksSENI[,6]==1
LS6[is.na(LS6)]<-FALSE

LS7<-linksSENI[,7]==1
LS7[is.na(LS7)]<-FALSE
linksSENI5[LS6,]<-0
linksSENI5[LS7,]<-NA

# environmental group engagement

linksENVI<-cbind(		DATA$Links.Environmental.Member,
						DATA$Links.Environmental.Activities,
						DATA$Links.Environmental.Donated,
						DATA$Links.Environmental.VolunteerWork,
						DATA$Links.Environmental.PaidWork,
						DATA$Links.Environmental.None,
						DATA$Links.Environmental.Don.t.know)

linksENVI5<-linksENVI[,1:5]
LE6<-linksENVI[,6]==1

# sixth column = "none"
# seventh column = "don't know"

LE6[is.na(LE6)]<-FALSE
LE7<-linksENVI[,7]==1
LE7[is.na(LE7)]<-FALSE
linksENVI5[LE6,]<-0
linksENVI5[LE7,]<-NA

# the order here is: ASYL (HURI), CCTV (HURI), CULT (CULT), CARE (SENI) and ROAD (ENVI)

rbound<-rbind(linksHURI5,linksHURI5,linksCULT5,linksSENI5,linksENVI5)

SL<-rowMeans(rbound,na.rm=TRUE)

SB<-1-as.numeric(SL==0)
SB[is.na(SL)]<-NA

# 2.3.2   Contacts
# 2.3.2.1 Business group contact

BusinessContacts<-DATA$Contacts.BusinessGroups
FirmContacts<-DATA$Contacts.Business

BusiFirm<-cbind(BusinessContacts,FirmContacts)

# footnote below Table A3
coefB<-coefH(BusiFirm[rowSums(is.na(BusiFirm))==0,])
coefB$H

BusiFirms<-rowMeans(BusiFirm)

BC<-c(BusiFirms,BusiFirms,BusiFirms,BusiFirms,BusiFirms)

# 2.3.2.2 Issue-specific, societal contacts

LinkASYL<-DATA$Contacts.HumanRights
LinkCCTV<-DATA$Contacts.HumanRights
LinkCULT<-DATA$Contacts.Culture
LinkCARE<-DATA$Contacts.Seniors
LinkROAD<-DATA$Contacts.Environmental

CSp<-c(LinkASYL,LinkCCTV,LinkCULT,LinkCARE,LinkROAD)

# 2.3.2.3 Societal group contacts (broad)

SocLinks<-cbind(DATA$Contacts.Environmental,DATA$Contacts.Seniors,DATA$Contacts.HumanRights,DATA$Contacts.Sport,DATA$Contacts.Culture,DATA$Contacts.Charity,DATA$Contacts.Consumer,DATA$Contacts.Minority,DATA$Contacts.Identity)

#below Table A3
coefS<-coefH(SocLinks[rowSums(is.na(SocLinks))==0,])
coefS$H

SocLinks2<-rowMeans(SocLinks)

CS1<-c(SocLinks2,SocLinks2,SocLinks2,SocLinks2,SocLinks2)

# 2.3.2.4 Diffuse, societal group contacts

DifLinks<-cbind(DATA$Contacts.Environmental,DATA$Contacts.Consumer,DATA$Contacts.HumanRights)

#below Table A3
coefD<-coefH(DifLinks[rowSums(is.na(DifLinks))==0,])
coefD$H

DifLinks2<-rowMeans(DifLinks)

CS2<-c(DifLinks2,DifLinks2,DifLinks2,DifLinks2,DifLinks2)

# 2.3.3   Information
# 2.3.3.1 Business information

BusinessInfo<-DATA$Information.BusinessGroups

BI<-c(BusinessInfo,BusinessInfo,BusinessInfo,BusinessInfo,BusinessInfo)

# 2.3.3.2 NGO information

NI<-c(DATA$Information.NGO,DATA$Information.NGO,DATA$Information.NGO,DATA$Information.NGO,DATA$Information.NGO)

# 2.3.3.3 Business Information –NGO Information

BusiNGO<- DATA$Information.BusinessGroups - DATA$Information.NGO

BN<-c(BusiNGO,BusiNGO,BusiNGO,BusiNGO,BusiNGO)

########################################################################################
# PART 2.4: EGOCENTRIC MEASURES
# There are six measures of egocentrism: extremism, congruence (weighted and unweighted), dichotomous congruence (weighted and unweighted) and own position

# 2.4.1   Extremism 

ExtASYL<-abs(DATA$Position.Asylumseekers-3)
ExtCCTV<-abs(DATA$Position.CCTV-3)
ExtCULT<-abs(DATA$Position.Culture-3)
ExtCARE<-abs(DATA$Position.ElderlyCare-3)
ExtROAD<-abs(DATA$Position.TrafficJams-3)
Ext<-c(ExtASYL,ExtCCTV,ExtCULT,ExtCARE,ExtROAD)

# 2.4.2   Congruence
# Congruence for weighted data

# these are identify cases with specific values, those vectors can't handle missing data
Sel3ASYL<-(DATA$Position.Asylumseekers>3)
Sel3CCTV<-(DATA$Position.CCTV>3)
Sel3CULT<-(DATA$Position.Culture>3)
Sel3CARE<-(DATA$Position.ElderlyCare>3)
Sel3ROAD<-(DATA$Position.TrafficJams>3)

Sel3ASYL[is.na(Sel3ASYL)]<-FALSE
Sel3CCTV[is.na(Sel3CCTV)]<-FALSE
Sel3CULT[is.na(Sel3CULT)]<-FALSE
Sel3CARE[is.na(Sel3CARE)]<-FALSE
Sel3ROAD[is.na(Sel3ROAD)]<-FALSE

Sel3<-c(Sel3ASYL,Sel3CCTV,Sel3CULT,Sel3CARE,Sel3ROAD)

Sel2ASYL<-(DATA$Position.Asylumseekers<4)
Sel2CCTV<-(DATA$Position.CCTV<4)
Sel2CULT<-(DATA$Position.Culture<4)
Sel2CARE<-(DATA$Position.ElderlyCare<4)
Sel2ROAD<-(DATA$Position.TrafficJams<4)
Sel2ASYL[is.na(Sel2ASYL)]<-FALSE
Sel2CCTV[is.na(Sel2CCTV)]<-FALSE
Sel2CULT[is.na(Sel2CULT)]<-FALSE
Sel2CARE[is.na(Sel2CARE)]<-FALSE
Sel2ROAD[is.na(Sel2ROAD)]<-FALSE

Sel2<-c(Sel2ASYL,Sel2CCTV,Sel2CULT,Sel2CARE,Sel2ROAD)

# these are the variables we are working with
ConASYL<-rep(NA,dim(DATA)[1])
ConCCTV<-rep(NA,dim(DATA)[1])
ConCULT<-rep(NA,dim(DATA)[1])
ConCARE<-rep(NA,dim(DATA)[1])
ConROAD<-rep(NA,dim(DATA)[1])

#assign values if answer is greater than 3, assign reverse if not
ConASYL[Sel3ASYL]<-TruASYL[Sel3ASYL]
ConASYL[Sel2ASYL]<-100-TruASYL[Sel2ASYL]

ConCCTV[Sel3CCTV]<-TruCCTV[Sel3CCTV]
ConCCTV[Sel2CCTV]<-100-TruCCTV[Sel2CCTV]

ConCULT[Sel3CULT]<-TruCULT[Sel3CULT]
ConCULT[Sel2CULT]<-100-TruCULT[Sel2CULT]

ConCARE[Sel3CARE]<-TruCARE[Sel3CARE]
ConCARE[Sel2CARE]<-100-TruCARE[Sel2CARE]

ConROAD[Sel3ROAD]<-TruROAD[Sel3ROAD]
ConROAD[Sel2ROAD]<-100-TruROAD[Sel2ROAD]

Con<-c(ConASYL,ConCCTV,ConCULT,ConCARE,ConROAD)

#divide by 100 for interpretation in some regressions
Con1<-Con/100

# Congruence for unweighted data

CwnASYL<-rep(NA,dim(DATA)[1])
CwnCCTV<-rep(NA,dim(DATA)[1])
CwnCULT<-rep(NA,dim(DATA)[1])
CwnCARE<-rep(NA,dim(DATA)[1])
CwnROAD<-rep(NA,dim(DATA)[1])

#assign values if answer is greater than 3, assign reverse if not
CwnASYL[Sel3ASYL]<-TrwASYL[Sel3ASYL]
CwnASYL[Sel2ASYL]<-100-TrwASYL[Sel2ASYL]

CwnCCTV[Sel3CCTV]<-TrwCCTV[Sel3CCTV]
CwnCCTV[Sel2CCTV]<-100-TrwCCTV[Sel2CCTV]

CwnCULT[Sel3CULT]<-TrwCULT[Sel3CULT]
CwnCULT[Sel2CULT]<-100-TrwCULT[Sel2CULT]

CwnCARE[Sel3CARE]<-TrwCARE[Sel3CARE]
CwnCARE[Sel2CARE]<-100-TrwCARE[Sel2CARE]

CwnROAD[Sel3ROAD]<-TrwROAD[Sel3ROAD]
CwnROAD[Sel2ROAD]<-100-TrwROAD[Sel2ROAD]

Cwn<-c(CwnASYL,CwnCCTV,CwnCULT,CwnCARE,CwnROAD)

#divide by 100 for interpretation in some regressions
Cwn1<-Cwn/100

# 2.4.3   Own position

#reverse for flipped positions
SelfLR<-c(6-DATA$Position.Asylumseekers,DATA$Position.CCTV,DATA$Position.Culture,6-DATA$Position.ElderlyCare,DATA$Position.TrafficJams)

########################################################################################
# PART 2.5: CONTROL VARIABLES

# 2.5.1   Expertise

SpcASYL<-DATA$Issue.IntegrationandImmigration
SpcCCTV<-DATA$Issue.JusticeandSafety
SpcCULT<-DATA$Issue.Culture
SpcCARE<-DATA$Issue.Healthcare
SpcROAD<-DATA$Issue.Transport

Spc<-c(SpcASYL,SpcCCTV,SpcCULT,SpcCARE,SpcROAD)

# 2.5.2   PPG Size

Seats<-DATA$NUM

Seat2<-c(Seats,Seats,Seats,Seats,Seats)

# 2.5.3   Levels

# Regional
LvlREG<-DATA$Level=="REG"

# National
LvlReg<-c(LvlREG,LvlREG,LvlREG,LvlREG,LvlREG)

LvlNAT<-DATA$Level=="NAT"

LvlNat<-c(LvlNAT,LvlNAT,LvlNAT,LvlNAT,LvlNAT)

# Country
Countr<-DATA$Country=="DK"

Country<-c(Countr,Countr,Countr,Countr,Countr)

# 2.5.4   Gender

FEMA<-DATA$Gender=="Female"

FEMAL<-c(FEMA,FEMA,FEMA,FEMA,FEMA)

# 2.5.5   Age

Age<-DATA$Background.Age

AGEL<-c(Age,Age,Age,Age,Age)

# 2.5.5   Education

BAMA<-DATA$Background.EducationBAMA

BAMAL<-c(BAMA,BAMA,BAMA,BAMA,BAMA)

# 2.5.5   Left-right position

LIRE<-DATA$Position.LeftRight

LIREL<-c(LIRE,LIRE,LIRE,LIRE,LIRE)

# 2.5.5   IDs for multilevel

N<-dim(DATA)[1]
ID1<-rep(1:N,5)
ID2<-c(rep(1,N),rep(2,N),rep(3,N),rep(4,N),rep(5,N))+Country*10

# 2.5.6   Priority

PriASYL<-DATA$Importance.Asylumseekers
PriCCTV<-DATA$Importance.CCTV
PriCULT<-DATA$Importance.Culture
PriCARE<-DATA$Importance.ElderlyCare
PriROAD<-DATA$Importance.TrafficJams 

Pri<-c(PriASYL,PriCCTV,PriCULT,PriCARE,PriROAD)

# 2.5.7	  Experience

EXPMC<-round(DATA$General.MunicipalCouncillor,0)
EXPRC<-round(DATA$General.RegionalCouncillor,0)
EXPNC<-round(DATA$General.NationalMP,0)

# missing value means no experience
EXPMC[is.na(EXPMC)]<-0
EXPRC[is.na(EXPRC)]<-0
EXPNC[is.na(EXPNC)]<-0

EXP<-EXPMC+EXPRC+EXPNC

EXP5<-c(EXP,EXP,EXP,EXP,EXP)

########################################################################################
# PART 3: DESCRIPTIVES

########################################################################################
# PART 3.1: TABLE 1

round(mean(DATA$Importance.Asylumseekers[DATA$Country=="DK"],na.rm=TRUE),2)
round(mean(DATA$Importance.CCTV[DATA$Country=="DK"],na.rm=TRUE),2)
round(mean(DATA$Importance.ElderlyCare[DATA$Country=="DK"],na.rm=TRUE),2)
round(mean(DATA$Importance.Culture[DATA$Country=="DK"],na.rm=TRUE),2)
round(mean(DATA$Importance.TrafficJams[DATA$Country=="DK"],na.rm=TRUE),2)

round(mean(DATA$Importance.Asylumseekers[DATA$Country=="NL"],na.rm=TRUE),2)
round(mean(DATA$Importance.CCTV[DATA$Country=="NL"],na.rm=TRUE),2)
round(mean(DATA$Importance.ElderlyCare[DATA$Country=="NL"],na.rm=TRUE),2)
round(mean(DATA$Importance.Culture[DATA$Country=="NL"],na.rm=TRUE),2)
round(mean(DATA$Importance.TrafficJams[DATA$Country=="NL"],na.rm=TRUE),2)

########################################################################################
# PART 3.2: TABLE 2

round(TruASYL[DATA$Country=="DK"][1],0)
round(TruASYL[DATA$Country=="NL"][1],0)
round(mean(EstASYL[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(EstASYL[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(AbsASYL[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(AbsASYL[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(RightASYL[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(RightASYL[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(ConASYL[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(ConASYL[DATA$Country=="NL"],na.rm=TRUE),0)

round(TruCCTV[DATA$Country=="DK"][1],0)
round(TruCCTV[DATA$Country=="NL"][1],0)
round(mean(EstCCTV[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(EstCCTV[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(AbsCCTV[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(AbsCCTV[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(RightCCTV[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(RightCCTV[DATA$Country=="NL"],na.rm=TRUE),0)

# to check values below 1: it should be written as <1% in the table
round(mean(RightCCTV[DATA$Country=="NL"],na.rm=TRUE),1)

round(mean(ConCCTV[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(ConCCTV[DATA$Country=="NL"],na.rm=TRUE),0)

round(TruCULT[DATA$Country=="DK"][1],0)
round(TruCULT[DATA$Country=="NL"][1],0)
round(mean(EstCULT[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(EstCULT[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(AbsCULT[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(AbsCULT[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(RightCULT[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(RightCULT[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(ConCULT[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(ConCULT[DATA$Country=="NL"],na.rm=TRUE),0)

round(TruCARE[DATA$Country=="DK"][1],0)
round(TruCARE[DATA$Country=="NL"][1],0)
round(mean(EstCARE[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(EstCARE[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(AbsCARE[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(AbsCARE[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(RightCARE[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(RightCARE[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(ConCARE[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(ConCARE[DATA$Country=="NL"],na.rm=TRUE),0)

round(TruROAD[DATA$Country=="DK"][1],0)
round(TruROAD[DATA$Country=="NL"][1],0)
round(mean(EstROAD[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(EstROAD[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(AbsROAD[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(AbsROAD[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(RightROAD[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(RightROAD[DATA$Country=="NL"],na.rm=TRUE),0)
round(mean(ConROAD[DATA$Country=="DK"],na.rm=TRUE),0)
round(mean(ConROAD[DATA$Country=="NL"],na.rm=TRUE),0)

# In text on pp.27-29
round(mean(AbsASYL,na.rm=TRUE),0)
round(mean(AbsCCTV,na.rm=TRUE),0)
round(mean(AbsCULT,na.rm=TRUE),0)
round(mean(AbsCARE,na.rm=TRUE),0)
round(mean(AbsROAD,na.rm=TRUE),0)

round(mean(Abs,na.rm=TRUE),0)
round(100-mean(Abs,na.rm=TRUE),0)

round(median(Abs,na.rm=TRUE),0)
round(mean(RightDif,na.rm=TRUE),0)

##########
# TABLE A58: 

round(mean(Abs[Country],na.rm=TRUE),0)
round(mean(Abs[!Country],na.rm=TRUE),0)

round(mean(RightDif[Country],na.rm=TRUE),0)
round(mean(RightDif[!Country],na.rm=TRUE),0)

round(mean(Dic[Country],na.rm=TRUE),2)
round(mean(Dic[!Country],na.rm=TRUE),2)

sum(Countr)
sum(!Countr)

########################################################################################
# PART 3.3: FIGURE 1

par( mfrow = c( 2, 3 ) )

hist(DifASYL,main="Histogram of Difference between Estimated and True
Share of Support on Asylum Question",xlab="Difference between Estimated and True
Share of Support on Asylum Question",ylim=c(0,475))

lines(x=c(0,0),y=c(0,1000),lwd=5)

hist(DifCCTV,main="Histogram of Difference between Estimated and True
Share of Support on Safety Question",xlab="Difference between Estimated and True
Share of Support on Safety Question",ylim=c(0,475))

lines(x=c(0,0),y=c(0,1000),lwd=5)

hist(DifCULT,main="Histogram of Difference between Estimated and True
Share of Support on Culture Question",xlab="Difference between Estimated and True
Share of Support on Culture Question",ylim=c(0,475))

lines(x=c(0,0),y=c(0,1000),lwd=5)

hist(DifCARE,main="Histogram of Difference between Estimated and True
Share of Support on Healthcare Question",xlab="Difference between Estimated and True
Share of Support on Healthcare Question",ylim=c(0,475))

lines(x=c(0,0),y=c(0,1000),lwd=5)

hist(DifROAD,main="Histogram of Difference between Estimated and True
Share of Support on Road Question",xlab="Difference between Estimated and True
Share of Support on Road Question",ylim=c(0,475))

lines(x=c(0,0),y=c(0,1000),lwd=5)

########################################################################################
# PART 3.4: TABLE A3

overview(Abs)
overview(RightDif)
overview(Dic)
overview(Ext)
overview(EXP5)
overview(Con)
overview(SelfLR)
overview(BB)
overview(SB)
overview(BC)
overview(CSp)
overview(CS1)
overview(CS2)
overview(NI)
overview(BI)
overview(BN)
overview(Spc)
overview(Seat2)
overview(EXP5)
overview(FEMAL)
overview(AGEL)
overview(BAMAL)
overview(LvlReg)
overview(LvlNat)
overview(Country)

########################################################################################
# PART 3.5: TABLE A4

cor.test(SB,EXP5)
cor.test(BB,EXP5)
cor.test(BC,EXP5)
cor.test(CSp,EXP5)
cor.test(CS1,EXP5)
cor.test(CS2,EXP5)
cor.test(BI,EXP5)
cor.test(NI,EXP5)
cor.test(BN,EXP5)
cor.test(SelfLR,EXP5)
cor.test(Con,EXP5)
cor.test(Ext,EXP5)

cor.test(Ext,BB)
cor.test(Ext,SB)
cor.test(Ext,BC)
cor.test(Ext,CSp)
cor.test(Ext,CS1)
cor.test(Ext,CS2)
cor.test(Ext,NI)
cor.test(Ext,BI)
cor.test(Ext,BN)
cor.test(Ext,SelfLR)
cor.test(Ext,Con)

cor.test(Con,BB)
cor.test(Con,SB)
cor.test(Con,BC)
cor.test(Con,CSp)
cor.test(Con,CS1)
cor.test(Con,CS2)
cor.test(Con,NI)
cor.test(Con,BI)
cor.test(Con,BN)
cor.test(Con,SelfLR)

cor.test(SelfLR,BB)
cor.test(SelfLR,SB)
cor.test(SelfLR,BC)
cor.test(SelfLR,CSp)
cor.test(SelfLR,CS1)
cor.test(SelfLR,CS2)
cor.test(SelfLR,NI)
cor.test(SelfLR,BI)
cor.test(SelfLR,BN)

cor.test(BN,BB)
cor.test(BN,SB)
cor.test(BN,BC)
cor.test(BN,CSp)
cor.test(BN,CS1)
cor.test(BN,CS2)
cor.test(BN,NI)
cor.test(BN,BI)

cor.test(NI,BB)
cor.test(NI,SB)
cor.test(NI,BC)
cor.test(NI,CSp)
cor.test(NI,CS1)
cor.test(NI,CS2)
cor.test(NI,BI)

cor.test(BI,BB)
cor.test(BI,SB)
cor.test(BI,BC)
cor.test(BI,CSp)
cor.test(BI,CS1)
cor.test(BI,CS2)

cor.test(CS2,BB)
cor.test(CS2,SB)
cor.test(CS2,BC)
cor.test(CS2,CSp)
cor.test(CS2,CS1)

cor.test(CS1,BB)
cor.test(CS1,SB)
cor.test(CS1,BC)
cor.test(CS1,CSp)

cor.test(CSp,BB)
cor.test(CSp,SB)
cor.test(CSp,BC)

cor.test(BC,BB)
cor.test(BC,SB)

cor.test(BB,SB)

########################################################################################
# PART 3.6: APPENDIX 1

DKEX<-(DATA$Minister=="1. VICEBORGMESTER")+(DATA$Minister=="2. VICEBORGMESTER")+(DATA$Minister=="BORGMESTER")+(DATA$Minister=="BØRNE- OG UNGDOMSBORGMESTER")+(DATA$Minister=="VICEBORGMESTER")

# as percentage
round(sum(DKEX,na.rm=TRUE)/sum(DATA$Country=="DK",na.rm=TRUE)*100,1)

########################################################################################
# PART 4: ANALYSES

# TABLE 3

model0001<-lmer(Abs~Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0002<-lmer(Abs~Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0003<-lmer(Abs~Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0004<-lmer(Abs~Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer001<-htmlreg(list(model0001,model0002,model0003,model0004))

write(stargazer001,"results/Table 3.html")

# numbers on p.30

# the effect of congruence
round(max(Con,na.rm=TRUE)*summary(model0003)[[10]][[2]]-min(Con,na.rm=TRUE)*summary(model0003)[[10]][[2]],0)

# the range of the effect of business groups between -3 and -4
round(max(BC,na.rm=TRUE)*summary(model0002)[[10]][[3]]-min(BC,na.rm=TRUE)*summary(model0002)[[10]][[3]],0)

round(max(BC,na.rm=TRUE)*summary(model0003)[[10]][[3]]-min(BC,na.rm=TRUE)*summary(model0003)[[10]][[3]],0)

round(max(BC,na.rm=TRUE)*summary(model0004)[[10]][[3]]-min(BC,na.rm=TRUE)*summary(model0004)[[10]][[3]],0)


# FIGURE A1
# change variable labels because these will be the labels in figure
Accuracy<-Abs
Congruence<-Con
Business <-BC
Societal <-CS1
Gender <- FEMAL
Age <- AGEL
Education <- BAMAL
Seats <- Seat2
Specialization <- Spc
Regional <- LvlReg
National <- LvlNat
Own <- SelfLR
Bias <- RightDif
Respondents<-ID1
Country_Issue<-ID2

# regenerate regression
model0003b<-lmer(Accuracy~Congruence+Business+Societal+Gender+Age+Education+Seats+Specialization+Regional+National+Country+(1|Respondents)+(1|Country_Issue))

# Figure A1 in order from top left then top right, to bottom right

# Figure A1 - 1
check_model(model0003b,check="pp_check")

# Figure A1 - 2
check_model(model0003b,check="linearity")

# Figure A1 - 3
check_model(model0003b,check="homogeneity")

# Figure A1 - 4
check_model(model0003b,check="outliers")

# Figure A1 - 5
check_model(model0003b,check="vif")

# Figure A1 - 6
check_model(model0003b,check="normality")

# Figure A1 - 7-8
check_model(model0003b,check="reqq")


# TABLE 4

model0101<-lmer(RightDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# This model reports a "boundary (singular) fit: see help('isSingular')". This indicates that the random effects variance-covariance matrix is of less than full rank. Upon investigation, we found that the random effects for respondent have zero variance. Bolker (2024) GLMM FAQ (https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#singular-models-random-effect-variances-estimated-as-zero-or-correlations-estimated-as---1) proposes a number of solutions. One is that "[i]f a variance component is zero, dropping it from the model will have no effect on any of the estimated quantities (although it will affect the AIC, as the variance parameter is counted even though it has no effect).” leading to the conclusion that “if one chooses for philosophical grounds to retain these parameters, it won’t change any of the answers." We ran an alternative model in which we dropped the random effects for respondents and compared the estimates for the coefficients. We found that the relevant coefficients and standard errors were identical to three decimal places (which is one more than reported in the paper). In the interest of consistency between the models predicting bias and accuracy, we decided to maintain the same  structure in terms of the random intercepts in the bias models as in the other models. Note that we do not focus on measures of fit such as the AIC in the paper .

model0101b<-lmer(RightDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

round((summary(model0101))$coefficients[,1]-(summary(model0101b))$coefficients[,1],3)
round((summary(model0101))$coefficients[,2]-(summary(model0101b))$coefficients[,2],3)

model0102<-lmer(RightDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# This model reports a "boundary (singular) fit: see help('isSingular')". As above we run the model without respondent random effects and found that the results are identical to three decimal places.

model0102b<-lmer(RightDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

round((summary(model0102))$coefficients[,1]-(summary(model0102b))$coefficients[,1],3)
round((summary(model0102))$coefficients[,2]-(summary(model0102b))$coefficients[,2],3)

model0103<-lmer(RightDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# This model reports a "boundary (singular) fit: see help('isSingular')". As above we run the model without respondent random effects and found that the results are identical to three decimal places.

model0103b<-lmer(RightDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

round((summary(model0103))$coefficients[,1]-(summary(model0103b))$coefficients[,1],3)
round((summary(model0103))$coefficients[,2]-(summary(model0103b))$coefficients[,2],3)

model0104<-lmer(RightDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# This model reports a "boundary (singular) fit: see help('isSingular')". As above we run the model without respondent random effects and found that the results are identical to three decimal places.

model0104b<-lmer(RightDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

round((summary(model0104))$coefficients[,1]-(summary(model0104b))$coefficients[,1],3)
round((summary(model0104))$coefficients[,2]-(summary(model0104b))$coefficients[,2],3)


stargazer004<-htmlreg(list(model0101,model0102,model0103,model0104))

write(stargazer004,"results/Table 4.html")

# numbers on p.35

# effect of own position
round(max(SelfLR,na.rm=TRUE)*summary(model0103)[[10]][[2]]-min(SelfLR,na.rm=TRUE)*summary(model0103)[[10]][[2]],0)

# the range of the effect of business groups is zero
round(max(BC,na.rm=TRUE)*summary(model0102)[[10]][[3]]-min(BC,na.rm=TRUE)*summary(model0102)[[10]][[3]],0)

round(max(BC,na.rm=TRUE)*summary(model0103)[[10]][[3]]-min(BC,na.rm=TRUE)*summary(model0103)[[10]][[3]],0)

round(max(BC,na.rm=TRUE)*summary(model0104)[[10]][[3]]-min(BC,na.rm=TRUE)*summary(model0104)[[10]][[3]],0)

# FIGURE A2
# regenerate regression

model0103b<-lmer(Bias~Own+Business+Societal+Gender+Age+Education+Seats+Specialization+Regional+National+Country+(1|Respondents)+(1|Country_Issue))

# Figure A2 in order from top left then top right, to bottom right

# Figure A2 - 1
check_model(model0103b,check="pp_check")

# Figure A1 - 2
check_model(model0103b,check="linearity")

# Figure A1 - 3
check_model(model0103b,check="homogeneity")

# Figure A1 - 4
check_model(model0103b,check="outliers")

# Figure A1 - 5
check_model(model0103b,check="vif")

# Figure A1 - 6
check_model(model0103b,check="normality")

# Figure A1 - 7-8
check_model(model0103b,check="reqq")

# TABLE A5

model0005<-lmer(Abs~Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0006<-lmer(Abs~Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0007<-lmer(Abs~Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0008<-lmer(Abs~Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer002<-htmlreg(list(model0005,model0006,model0007,model0008))

write(stargazer002,"results/Table A05.html")

# TABLE A6

model0010<-lmer(Abs~Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0011<-lmer(Abs~Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0012<-lmer(Abs~Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0013<-lmer(Abs~Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0014<-lmer(Abs~Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0015<-lmer(Abs~Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0016<-lmer(Abs~Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer003<-htmlreg(list(model0010,model0011,model0012,model0013,model0014,model0015,model0016))

write(stargazer003,"results/Table A06.html")

# TABLE A7

model0201<-lmer(Abs~Ext+Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0202<-lmer(Abs~Ext+Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0203<-lmer(Abs~Ext+Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0204<-lmer(Abs~Ext+Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0205<-lmer(Abs~Ext+Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0206<-lmer(Abs~Ext+Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0207<-lmer(Abs~Ext+Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0208<-lmer(Abs~Ext+Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer007<-htmlreg(list(model0201,model0202,model0203,model0204,model0205,model0206,model0207,model0208))

write(stargazer007,"results/Table A07.html")

# TABLE A8

model0209<-lmer(Abs~Ext+Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0210<-lmer(Abs~Ext+Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0211<-lmer(Abs~Ext+Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0212<-lmer(Abs~Ext+Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0213<-lmer(Abs~Ext+Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0214<-lmer(Abs~Ext+Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0215<-lmer(Abs~Ext+Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer008<-htmlreg(list(model0209,model0210,model0211,model0212,model0213,model0214,model0215))

write(stargazer008,"results/Table A08.html")

# TABLE A9

model0301<-lmer(Abs~FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0302<-lmer(Abs~BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0303<-lmer(Abs~BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0304<-lmer(Abs~BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0305<-lmer(Abs~BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0306<-lmer(Abs~CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0307<-lmer(Abs~CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0308<-lmer(Abs~CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer009<-htmlreg(list(model0301,model0302,model0303,model0304,model0305,model0306,model0307,model0308))

write(stargazer009,"results/Table A09.html")

# TABLE A10

model0309<-lmer(Abs~BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0310<-lmer(Abs~SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0311<-lmer(Abs~BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0312<-lmer(Abs~BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0313<-lmer(Abs~NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0314<-lmer(Abs~BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0315<-lmer(Abs~BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer010<-htmlreg(list(model0309,model0310,model0311,model0312,model0313,model0314,model0315))

write(stargazer010,"results/Table A10.html")

#####

# TABLE A11

model0401<-lmer(Aws~Cwn+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0402<-lmer(Aws~Cwn+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0403<-lmer(Aws~Cwn+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0404<-lmer(Aws~Cwn+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0405<-lmer(Aws~Cwn+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0406<-lmer(Aws~Cwn+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0407<-lmer(Aws~Cwn+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0408<-lmer(Aws~Cwn+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer019<-htmlreg(list(model0401,model0402,model0403,model0404,model0405,model0406,model0407,model0408))

write(stargazer019,"results/Table A11.html")

# TABLE A12

model0409<-lmer(Aws~Cwn+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0410<-lmer(Aws~Cwn+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0411<-lmer(Aws~Cwn+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0412<-lmer(Aws~Cwn+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0413<-lmer(Aws~Cwn+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0414<-lmer(Aws~Cwn+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0415<-lmer(Aws~Cwn+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer020<-htmlreg(list(model0409,model0410,model0411,model0412,model0413,model0414,model0415))

write(stargazer020,"results/Table A12.html")

#####

# TABLE A13

model0701<-lmer(Abs~Pri+Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0702<-lmer(Abs~Pri+Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0703<-lmer(Abs~Pri+Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0704<-lmer(Abs~Pri+Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0705<-lmer(Abs~Pri+Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0706<-lmer(Abs~Pri+Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0707<-lmer(Abs~Pri+Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0708<-lmer(Abs~Pri+Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer107<-htmlreg(list(model0701,model0702,model0703,model0704,model0705,model0706,model0707,model0708))

write(stargazer107,"results/Table A13.html")

# TABLE A14

model0709<-lmer(Abs~Pri+Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0710<-lmer(Abs~Pri+Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0711<-lmer(Abs~Pri+Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0712<-lmer(Abs~Pri+Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0713<-lmer(Abs~Pri+Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0714<-lmer(Abs~Pri+Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0715<-lmer(Abs~Pri+Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer108<-htmlreg(list(model0709,model0710,model0711,model0712,model0713,model0714,model0715))

write(stargazer108,"results/Table A14.html")

#####
# TABLE A15

model0801<-lmer(Abs~EXP5+Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0802<-lmer(Abs~EXP5+Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0803<-lmer(Abs~EXP5+Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0804<-lmer(Abs~EXP5+Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0805<-lmer(Abs~EXP5+Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0806<-lmer(Abs~EXP5+Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0807<-lmer(Abs~EXP5+Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0808<-lmer(Abs~EXP5+Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer807<-htmlreg(list(model0801,model0802,model0803,model0804,model0805,model0806,model0807,model0808))

write(stargazer807,"results/Table A15.html")

# TABLE A16

model0809<-lmer(Abs~EXP5+Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0810<-lmer(Abs~EXP5+Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0811<-lmer(Abs~EXP5+Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0812<-lmer(Abs~EXP5+Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0813<-lmer(Abs~EXP5+Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0814<-lmer(Abs~EXP5+Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0815<-lmer(Abs~EXP5+Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer808<-htmlreg(list(model0809,model0810,model0811,model0812,model0813,model0814,model0815))

write(stargazer808,"results/Table A16.html")

# TABLE A17
# split the data by country
# country = 1 is Denmark; country = 0 is the Netherlands

NLDK<-as.data.frame(as.matrix(cbind(Abs,Con,FEMAL,AGEL,BAMAL,Seat2,Spc,LvlReg,LvlNat,Country,ID1,ID2,BB,SB,BC,CSp,CS1,CS2,BI,NI,BN)))

modelnl01<-lmer(Abs~Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl02<-lmer(Abs~Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl03<-lmer(Abs~Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl04<-lmer(Abs~Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl05<-lmer(Abs~Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl06<-lmer(Abs~Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl07<-lmer(Abs~Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl08<-lmer(Abs~Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])

stargazernl01<-htmlreg(list(modelnl01,modelnl02,modelnl03,modelnl04,modelnl05,modelnl06,modelnl07,modelnl08))

write(stargazernl01,"results/Table A17.html")

# TABLE A18

modelnl09<-lmer(Abs~Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl10<-lmer(Abs~Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl11<-lmer(Abs~Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl12<-lmer(Abs~Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl13<-lmer(Abs~Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl14<-lmer(Abs~Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl15<-lmer(Abs~Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])

stargazernl02<-htmlreg(list(modelnl09,modelnl10,modelnl11,modelnl12,modelnl13,modelnl14,modelnl15))

write(stargazernl02,"results/Table A18.html")

# TABLE A19

modelDK01<-lmer(Abs~Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK02<-lmer(Abs~Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK03<-lmer(Abs~Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK04<-lmer(Abs~Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK05<-lmer(Abs~Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK06<-lmer(Abs~Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK07<-lmer(Abs~Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK08<-lmer(Abs~Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])

stargazerDK01<-htmlreg(list(modelDK01,modelDK02,modelDK03,modelDK04,modelDK05,modelDK06,modelDK07,modelDK08))

write(stargazerDK01,"results/Table A19.html")

# TABLE A20

modelDK09<-lmer(Abs~Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK10<-lmer(Abs~Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK11<-lmer(Abs~Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK12<-lmer(Abs~Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK13<-lmer(Abs~Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK14<-lmer(Abs~Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modelDK15<-lmer(Abs~Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])

stargazerDK02<-htmlreg(list(modelDK09,modelDK10,modelDK11,modelDK12,modelDK13,modelDK14,modelDK15))

write(stargazerDK02,"results/Table A20.html")

# TABLE A21
# uses transformed accuracy

model9901<-lmer(Sba~Con+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9902<-lmer(Sba~Con+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9903<-lmer(Sba~Con+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9904<-lmer(Sba~Con+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9905<-lmer(Sba~Con+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9906<-lmer(Sba~Con+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9907<-lmer(Sba~Con+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9908<-lmer(Sba~Con+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer923<-htmlreg(list(model9901,model9902,model9903,model9904,model9905,model9906,model9907,model9908))

write(stargazer923,"results/Table A21.html")

# TABLE A22
# uses transformed accuracy

model9909<-lmer(Sba~Con+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9910<-lmer(Sba~Con+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9911<-lmer(Sba~Con+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9912<-lmer(Sba~Con+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9913<-lmer(Sba~Con+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9914<-lmer(Sba~Con+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model9915<-lmer(Sba~Con+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer924<-htmlreg(list(model9909,model9910,model9911,model9912,model9913,model9914,model9915))

write(stargazer924,"results/Table A22.html")

# FIGURE A3
# change variable labels because these will be the labels in figure

Trans.Accuracy<-Sba
Congruence<-Con
Business <-BC
Societal <-CS1
Gender <- FEMAL
Age <- AGEL
Education <- BAMAL
Seats <- Seat2
Specialization <- Spc
Regional <- LvlReg
National <- LvlNat
Own <- SelfLR
Bias <- RightDif

# regenerate regression
model9903b<-lmer(Trans.Accuracy~Congruence+Business+Societal+Gender+Age+Education+Seats+Specialization+Regional+National+Country+(1|Respondents)+(1|Country_Issue))

# Figure A3 in order from top left then top right, to bottom right

# Figure A3 - 1
check_model(model9903b,check="pp_check")

# Figure A3 - 2
check_model(model9903b,check="linearity")

# Figure A3 - 3
check_model(model9903b,check="homogeneity")

# Figure A3 - 4
check_model(model9903b,check="outliers")

# Figure A3 - 5
check_model(model9903b,check="vif")

# Figure A3 - 6
check_model(model9903b,check="normality")

# Figure A3 - 7-8
check_model(model9903b,check="reqq")

# TABLE A23
# in the beta regression we use the accuracy and congruence divided by one hundred 

model8801<-glmmTMB(Abs1~Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8802<-glmmTMB(Abs1~Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8803<-glmmTMB(Abs1~Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8804<-glmmTMB(Abs1~Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8805<-glmmTMB(Abs1~Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8806<-glmmTMB(Abs1~Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8807<-glmmTMB(Abs1~Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8808<-glmmTMB(Abs1~Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))

write(htmlreg(list(model8801,model8802,model8803,model8804,model8805,model8806,model8807,model8808)),"results/Table A23.html")

# TABLE A24

model8809<-glmmTMB(Abs1~Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8810<-glmmTMB(Abs1~Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8811<-glmmTMB(Abs1~Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8812<-glmmTMB(Abs1~Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8813<-glmmTMB(Abs1~Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8814<-glmmTMB(Abs1~Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))
model8815<-glmmTMB(Abs1~Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),family=beta_family(link = "logit"))

write(htmlreg(list(model8809,model8810,model8811,model8812,model8813,model8814,model8815)),"results/Table A24.html")

# TABLE A25

model0105<-lmer(RightDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0106<-lmer(RightDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0107<-lmer(RightDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0108<-lmer(RightDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

model0105b<-lmer(RightDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0106b<-lmer(RightDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0107b<-lmer(RightDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0108b<-lmer(RightDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.

round((summary(model0105))$coefficients[,1]-(summary(model0105b))$coefficients[,1],3)
round((summary(model0106))$coefficients[,1]-(summary(model0106b))$coefficients[,1],3)
round((summary(model0107))$coefficients[,1]-(summary(model0107b))$coefficients[,1],3)
round((summary(model0108))$coefficients[,1]-(summary(model0108b))$coefficients[,1],3)

round((summary(model0105))$coefficients[,2]-(summary(model0105b))$coefficients[,2],3)
round((summary(model0106))$coefficients[,2]-(summary(model0106b))$coefficients[,2],3)
round((summary(model0107))$coefficients[,2]-(summary(model0107b))$coefficients[,2],3)
round((summary(model0108))$coefficients[,2]-(summary(model0108b))$coefficients[,2],3)

stargazer005<-htmlreg(list(model0105,model0106,model0107,model0108))

write(stargazer005,"results/Table A25.html")

# TABLE A26

model0109<-lmer(RightDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0110<-lmer(RightDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0111<-lmer(RightDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0112<-lmer(RightDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0113<-lmer(RightDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0114<-lmer(RightDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0115<-lmer(RightDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.

model0109b<-lmer(RightDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0110b<-lmer(RightDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0111b<-lmer(RightDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0112b<-lmer(RightDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0113b<-lmer(RightDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0114b<-lmer(RightDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0115b<-lmer(RightDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))


round((summary(model0109))$coefficients[,1]-(summary(model0109b))$coefficients[,1],3)
round((summary(model0110))$coefficients[,1]-(summary(model0110b))$coefficients[,1],3)
round((summary(model0111))$coefficients[,1]-(summary(model0111b))$coefficients[,1],3)
round((summary(model0112))$coefficients[,1]-(summary(model0112b))$coefficients[,1],3)
round((summary(model0113))$coefficients[,1]-(summary(model0113b))$coefficients[,1],3)
round((summary(model0114))$coefficients[,1]-(summary(model0114b))$coefficients[,1],3)
round((summary(model0115))$coefficients[,1]-(summary(model0115b))$coefficients[,1],3)

round((summary(model0109))$coefficients[,2]-(summary(model0109b))$coefficients[,2],3)
round((summary(model0110))$coefficients[,2]-(summary(model0110b))$coefficients[,2],3)
round((summary(model0111))$coefficients[,2]-(summary(model0111b))$coefficients[,2],3)
round((summary(model0112))$coefficients[,2]-(summary(model0112b))$coefficients[,2],3)
round((summary(model0113))$coefficients[,2]-(summary(model0113b))$coefficients[,2],3)
round((summary(model0114))$coefficients[,2]-(summary(model0114b))$coefficients[,2],3)
round((summary(model0115))$coefficients[,2]-(summary(model0115b))$coefficients[,2],3)


stargazer006<-htmlreg(list(model0109,model0110,model0111,model0112,model0113,model0114,model0115))

write(stargazer006,"results/Table A26.html")

# TABLE A27

model0701<-lmer(RightDif~FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0702<-lmer(RightDif~BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0703<-lmer(RightDif~BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0704<-lmer(RightDif~BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0705<-lmer(RightDif~BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0706<-lmer(RightDif~CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0707<-lmer(RightDif~CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0708<-lmer(RightDif~CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.

model0701b<-lmer(RightDif~FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0702b<-lmer(RightDif~BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0703b<-lmer(RightDif~BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0704b<-lmer(RightDif~BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0705b<-lmer(RightDif~BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0706b<-lmer(RightDif~CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0707b<-lmer(RightDif~CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0708b<-lmer(RightDif~CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

round((summary(model0701))$coefficients[,1]-(summary(model0701b))$coefficients[,1],3)
round((summary(model0702))$coefficients[,1]-(summary(model0702b))$coefficients[,1],3)
round((summary(model0703))$coefficients[,1]-(summary(model0703b))$coefficients[,1],3)
round((summary(model0704))$coefficients[,1]-(summary(model0704b))$coefficients[,1],3)
round((summary(model0705))$coefficients[,1]-(summary(model0705b))$coefficients[,1],3)
round((summary(model0706))$coefficients[,1]-(summary(model0706b))$coefficients[,1],3)
round((summary(model0707))$coefficients[,1]-(summary(model0707b))$coefficients[,1],3)
round((summary(model0708))$coefficients[,1]-(summary(model0708b))$coefficients[,1],3)

round((summary(model0701))$coefficients[,2]-(summary(model0701b))$coefficients[,2],3)
round((summary(model0702))$coefficients[,2]-(summary(model0702b))$coefficients[,2],3)
round((summary(model0703))$coefficients[,2]-(summary(model0703b))$coefficients[,2],3)
round((summary(model0704))$coefficients[,2]-(summary(model0704b))$coefficients[,2],3)
round((summary(model0705))$coefficients[,2]-(summary(model0705b))$coefficients[,2],3)
round((summary(model0706))$coefficients[,2]-(summary(model0706b))$coefficients[,2],3)
round((summary(model0707))$coefficients[,2]-(summary(model0707b))$coefficients[,2],3)
round((summary(model0708))$coefficients[,2]-(summary(model0708b))$coefficients[,2],3)

stargazer017<-htmlreg(list(model0701,model0702,model0703,model0704,model0705,model0706,model0707,model0708))

write(stargazer017,"results/Table A27.html")

# TABLE A28

model0709<-lmer(RightDif~BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0710<-lmer(RightDif~SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0711<-lmer(RightDif~BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0712<-lmer(RightDif~BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0713<-lmer(RightDif~NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0714<-lmer(RightDif~BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0715<-lmer(RightDif~BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

model0709b<-lmer(RightDif~BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0710b<-lmer(RightDif~SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0711b<-lmer(RightDif~BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0712b<-lmer(RightDif~BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0713b<-lmer(RightDif~NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0714b<-lmer(RightDif~BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0715b<-lmer(RightDif~BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


round((summary(model0709))$coefficients[,1]-(summary(model0709b))$coefficients[,1],3)
round((summary(model0710))$coefficients[,1]-(summary(model0710b))$coefficients[,1],3)
round((summary(model0711))$coefficients[,1]-(summary(model0711b))$coefficients[,1],3)
round((summary(model0712))$coefficients[,1]-(summary(model0712b))$coefficients[,1],3)
round((summary(model0713))$coefficients[,1]-(summary(model0713b))$coefficients[,1],3)
round((summary(model0714))$coefficients[,1]-(summary(model0714b))$coefficients[,1],3)
round((summary(model0715))$coefficients[,1]-(summary(model0715b))$coefficients[,1],3)

round((summary(model0709))$coefficients[,2]-(summary(model0709b))$coefficients[,2],3)
round((summary(model0710))$coefficients[,2]-(summary(model0710b))$coefficients[,2],3)
round((summary(model0711))$coefficients[,2]-(summary(model0711b))$coefficients[,2],3)
round((summary(model0712))$coefficients[,2]-(summary(model0712b))$coefficients[,2],3)
round((summary(model0713))$coefficients[,2]-(summary(model0713b))$coefficients[,2],3)
round((summary(model0714))$coefficients[,2]-(summary(model0714b))$coefficients[,2],3)
round((summary(model0715))$coefficients[,2]-(summary(model0715b))$coefficients[,2],3)

stargazer018<-htmlreg(list(model0709,model0710,model0711,model0712,model0713,model0714,model0715))

write(stargazer018,"results/Table A28.html")

#####

# TABLE A29

model0901<-lmer(RwghtDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0902<-lmer(RwghtDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0903<-lmer(RwghtDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0904<-lmer(RwghtDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0905<-lmer(RwghtDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0906<-lmer(RwghtDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0907<-lmer(RwghtDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0908<-lmer(RwghtDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

model0901b<-lmer(RwghtDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0902b<-lmer(RwghtDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0903b<-lmer(RwghtDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0904b<-lmer(RwghtDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0905b<-lmer(RwghtDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0906b<-lmer(RwghtDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0907b<-lmer(RwghtDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0908b<-lmer(RwghtDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


round((summary(model0901))$coefficients[,1]-(summary(model0901b))$coefficients[,1],3)
round((summary(model0902))$coefficients[,1]-(summary(model0902b))$coefficients[,1],3)
round((summary(model0903))$coefficients[,1]-(summary(model0903b))$coefficients[,1],3)
round((summary(model0904))$coefficients[,1]-(summary(model0904b))$coefficients[,1],3)
round((summary(model0905))$coefficients[,1]-(summary(model0905b))$coefficients[,1],3)
round((summary(model0906))$coefficients[,1]-(summary(model0906b))$coefficients[,1],3)
round((summary(model0907))$coefficients[,1]-(summary(model0907b))$coefficients[,1],3)
round((summary(model0908))$coefficients[,1]-(summary(model0908b))$coefficients[,1],3)

round((summary(model0901))$coefficients[,2]-(summary(model0901b))$coefficients[,2],3)
round((summary(model0902))$coefficients[,2]-(summary(model0902b))$coefficients[,2],3)
round((summary(model0903))$coefficients[,2]-(summary(model0903b))$coefficients[,2],3)
round((summary(model0904))$coefficients[,2]-(summary(model0904b))$coefficients[,2],3)
round((summary(model0905))$coefficients[,2]-(summary(model0905b))$coefficients[,2],3)
round((summary(model0906))$coefficients[,2]-(summary(model0906b))$coefficients[,2],3)
round((summary(model0907))$coefficients[,2]-(summary(model0907b))$coefficients[,2],3)
round((summary(model0908))$coefficients[,2]-(summary(model0908b))$coefficients[,2],3)


stargazer021<-htmlreg(list(model0901,model0902,model0903,model0904,model0905,model0906,model0907,model0908))



write(stargazer021,"results/Table A29.html")

# TABLE A30

model0909<-lmer(RwghtDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0910<-lmer(RwghtDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0911<-lmer(RwghtDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0912<-lmer(RwghtDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0913<-lmer(RwghtDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0914<-lmer(RwghtDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0915<-lmer(RwghtDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

model0909b<-lmer(RwghtDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0910b<-lmer(RwghtDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0911b<-lmer(RwghtDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0912b<-lmer(RwghtDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0913b<-lmer(RwghtDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0914b<-lmer(RwghtDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model0915b<-lmer(RwghtDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.



round((summary(model0909))$coefficients[,1]-(summary(model0909b))$coefficients[,1],3)
round((summary(model0910))$coefficients[,1]-(summary(model0910b))$coefficients[,1],3)
round((summary(model0911))$coefficients[,1]-(summary(model0911b))$coefficients[,1],3)
round((summary(model0912))$coefficients[,1]-(summary(model0912b))$coefficients[,1],3)
round((summary(model0913))$coefficients[,1]-(summary(model0913b))$coefficients[,1],3)
round((summary(model0914))$coefficients[,1]-(summary(model0914b))$coefficients[,1],3)
round((summary(model0915))$coefficients[,1]-(summary(model0915b))$coefficients[,1],3)

round((summary(model0909))$coefficients[,2]-(summary(model0909b))$coefficients[,2],3)
round((summary(model0910))$coefficients[,2]-(summary(model0910b))$coefficients[,2],3)
round((summary(model0911))$coefficients[,2]-(summary(model0911b))$coefficients[,2],3)
round((summary(model0912))$coefficients[,2]-(summary(model0912b))$coefficients[,2],3)
round((summary(model0913))$coefficients[,2]-(summary(model0913b))$coefficients[,2],3)
round((summary(model0914))$coefficients[,2]-(summary(model0914b))$coefficients[,2],3)
round((summary(model0915))$coefficients[,2]-(summary(model0915b))$coefficients[,2],3)


stargazer022<-htmlreg(list(model0909,model0910,model0911,model0912,model0913,model0914,model0915))

write(stargazer022,"results/Table A30.html")

#####

# TABLE A31

model1301<-lmer(RightDif~SelfLR+Pri+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1302<-lmer(RightDif~SelfLR+Pri+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1303<-lmer(RightDif~SelfLR+Pri+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1304<-lmer(RightDif~SelfLR+Pri+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1305<-lmer(RightDif~SelfLR+Pri+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1306<-lmer(RightDif~SelfLR+Pri+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1307<-lmer(RightDif~SelfLR+Pri+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1308<-lmer(RightDif~SelfLR+Pri+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


model1301b<-lmer(RightDif~SelfLR+Pri+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1302b<-lmer(RightDif~SelfLR+Pri+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1303b<-lmer(RightDif~SelfLR+Pri+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1304b<-lmer(RightDif~SelfLR+Pri+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1305b<-lmer(RightDif~SelfLR+Pri+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1306b<-lmer(RightDif~SelfLR+Pri+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1307b<-lmer(RightDif~SelfLR+Pri+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1308b<-lmer(RightDif~SelfLR+Pri+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

round((summary(model1301))$coefficients[,1]-(summary(model1301b))$coefficients[,1],3)
round((summary(model1302))$coefficients[,1]-(summary(model1302b))$coefficients[,1],3)
round((summary(model1303))$coefficients[,1]-(summary(model1303b))$coefficients[,1],3)
round((summary(model1304))$coefficients[,1]-(summary(model1304b))$coefficients[,1],3)
round((summary(model1305))$coefficients[,1]-(summary(model1305b))$coefficients[,1],3)
round((summary(model1306))$coefficients[,1]-(summary(model1306b))$coefficients[,1],3)
round((summary(model1307))$coefficients[,1]-(summary(model1307b))$coefficients[,1],3)
round((summary(model1308))$coefficients[,1]-(summary(model1308b))$coefficients[,1],3)

round((summary(model1301))$coefficients[,2]-(summary(model1301b))$coefficients[,2],3)
round((summary(model1302))$coefficients[,2]-(summary(model1302b))$coefficients[,2],3)
round((summary(model1303))$coefficients[,2]-(summary(model1303b))$coefficients[,2],3)
round((summary(model1304))$coefficients[,2]-(summary(model1304b))$coefficients[,2],3)
round((summary(model1305))$coefficients[,2]-(summary(model1305b))$coefficients[,2],3)
round((summary(model1306))$coefficients[,2]-(summary(model1306b))$coefficients[,2],3)
round((summary(model1307))$coefficients[,2]-(summary(model1307b))$coefficients[,2],3)
round((summary(model1308))$coefficients[,2]-(summary(model1308b))$coefficients[,2],3)



stargazer029<-htmlreg(list(model1301,model1302,model1303,model1304,model1305,model1306,model1307,model1308))

write(stargazer029,"results/Table A31.html")

# TABLE A32

model1309<-lmer(RightDif~SelfLR+Pri+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1310<-lmer(RightDif~SelfLR+Pri+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1311<-lmer(RightDif~SelfLR+Pri+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1312<-lmer(RightDif~SelfLR+Pri+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1313<-lmer(RightDif~SelfLR+Pri+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1314<-lmer(RightDif~SelfLR+Pri+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1315<-lmer(RightDif~SelfLR+Pri+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

model1309b<-lmer(RightDif~SelfLR+Pri+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1310b<-lmer(RightDif~SelfLR+Pri+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1311b<-lmer(RightDif~SelfLR+Pri+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1312b<-lmer(RightDif~SelfLR+Pri+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1313b<-lmer(RightDif~SelfLR+Pri+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1314b<-lmer(RightDif~SelfLR+Pri+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1315b<-lmer(RightDif~SelfLR+Pri+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.



round((summary(model1309))$coefficients[,1]-(summary(model1309b))$coefficients[,1],3)
round((summary(model1310))$coefficients[,1]-(summary(model1310b))$coefficients[,1],3)
round((summary(model1311))$coefficients[,1]-(summary(model1311b))$coefficients[,1],3)
round((summary(model1312))$coefficients[,1]-(summary(model1312b))$coefficients[,1],3)
round((summary(model1313))$coefficients[,1]-(summary(model1313b))$coefficients[,1],3)
round((summary(model1314))$coefficients[,1]-(summary(model1314b))$coefficients[,1],3)
round((summary(model1315))$coefficients[,1]-(summary(model1315b))$coefficients[,1],3)

round((summary(model1309))$coefficients[,2]-(summary(model1309b))$coefficients[,2],3)
round((summary(model1310))$coefficients[,2]-(summary(model1310b))$coefficients[,2],3)
round((summary(model1311))$coefficients[,2]-(summary(model1311b))$coefficients[,2],3)
round((summary(model1312))$coefficients[,2]-(summary(model1312b))$coefficients[,2],3)
round((summary(model1313))$coefficients[,2]-(summary(model1313b))$coefficients[,2],3)
round((summary(model1314))$coefficients[,2]-(summary(model1314b))$coefficients[,2],3)
round((summary(model1315))$coefficients[,2]-(summary(model1315b))$coefficients[,2],3)

stargazer030<-htmlreg(list(model1309,model1310,model1311,model1312,model1313,model1314,model1315))

write(stargazer030,"results/Table A32.html")

######

# TABLE A33

model1401<-lmer(RightDif~SelfLR+EXP5+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1402<-lmer(RightDif~SelfLR+EXP5+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1403<-lmer(RightDif~SelfLR+EXP5+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1404<-lmer(RightDif~SelfLR+EXP5+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1405<-lmer(RightDif~SelfLR+EXP5+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1406<-lmer(RightDif~SelfLR+EXP5+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1407<-lmer(RightDif~SelfLR+EXP5+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
# to deal with non-convergence another optimizer was used
model1408<-lmer(RightDif~SelfLR+EXP5+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2),control = lmerControl(optimizer ="Nelder_Mead"))

model1401b<-lmer(RightDif~SelfLR+EXP5+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1402b<-lmer(RightDif~SelfLR+EXP5+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1403b<-lmer(RightDif~SelfLR+EXP5+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1404b<-lmer(RightDif~SelfLR+EXP5+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1405b<-lmer(RightDif~SelfLR+EXP5+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1406b<-lmer(RightDif~SelfLR+EXP5+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1407b<-lmer(RightDif~SelfLR+EXP5+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
# same change in optimizer as above
model1408b<-lmer(RightDif~SelfLR+EXP5+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2),control = lmerControl(optimizer ="Nelder_Mead"))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


round((summary(model1401))$coefficients[,1]-(summary(model1401b))$coefficients[,1],3)
round((summary(model1402))$coefficients[,1]-(summary(model1402b))$coefficients[,1],3)
round((summary(model1403))$coefficients[,1]-(summary(model1403b))$coefficients[,1],3)
round((summary(model1404))$coefficients[,1]-(summary(model1404b))$coefficients[,1],3)
round((summary(model1405))$coefficients[,1]-(summary(model1405b))$coefficients[,1],3)
round((summary(model1406))$coefficients[,1]-(summary(model1406b))$coefficients[,1],3)
round((summary(model1407))$coefficients[,1]-(summary(model1407b))$coefficients[,1],3)
round((summary(model1408))$coefficients[,1]-(summary(model1408b))$coefficients[,1],3)

round((summary(model1401))$coefficients[,2]-(summary(model1401b))$coefficients[,2],3)
round((summary(model1402))$coefficients[,2]-(summary(model1402b))$coefficients[,2],3)
round((summary(model1403))$coefficients[,2]-(summary(model1403b))$coefficients[,2],3)
round((summary(model1404))$coefficients[,2]-(summary(model1404b))$coefficients[,2],3)
round((summary(model1405))$coefficients[,2]-(summary(model1405b))$coefficients[,2],3)
round((summary(model1406))$coefficients[,2]-(summary(model1406b))$coefficients[,2],3)
round((summary(model1407))$coefficients[,2]-(summary(model1407b))$coefficients[,2],3)
round((summary(model1408))$coefficients[,2]-(summary(model1408b))$coefficients[,2],3)

stargazer033<-htmlreg(list(model1401,model1402,model1403,model1404,model1405,model1406,model1407,model1408))

write(stargazer033,"results/Table A33.html")

# TABLE A34

model1409<-lmer(RightDif~SelfLR+EXP5+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1410<-lmer(RightDif~SelfLR+EXP5+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1411<-lmer(RightDif~SelfLR+EXP5+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1412<-lmer(RightDif~SelfLR+EXP5+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1413<-lmer(RightDif~SelfLR+EXP5+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1414<-lmer(RightDif~SelfLR+EXP5+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1415<-lmer(RightDif~SelfLR+EXP5+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

model1409b<-lmer(RightDif~SelfLR+EXP5+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1410b<-lmer(RightDif~SelfLR+EXP5+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1411b<-lmer(RightDif~SelfLR+EXP5+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1412b<-lmer(RightDif~SelfLR+EXP5+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1413b<-lmer(RightDif~SelfLR+EXP5+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1414b<-lmer(RightDif~SelfLR+EXP5+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))
model1415b<-lmer(RightDif~SelfLR+EXP5+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID2))

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.

round((summary(model1409))$coefficients[,1]-(summary(model1409b))$coefficients[,1],3)
round((summary(model1410))$coefficients[,1]-(summary(model1410b))$coefficients[,1],3)
round((summary(model1411))$coefficients[,1]-(summary(model1411b))$coefficients[,1],3)
round((summary(model1412))$coefficients[,1]-(summary(model1412b))$coefficients[,1],3)
round((summary(model1413))$coefficients[,1]-(summary(model1413b))$coefficients[,1],3)
round((summary(model1414))$coefficients[,1]-(summary(model1414b))$coefficients[,1],3)
round((summary(model1415))$coefficients[,1]-(summary(model1415b))$coefficients[,1],3)

round((summary(model1409))$coefficients[,2]-(summary(model1409b))$coefficients[,2],3)
round((summary(model1410))$coefficients[,2]-(summary(model1410b))$coefficients[,2],3)
round((summary(model1411))$coefficients[,2]-(summary(model1411b))$coefficients[,2],3)
round((summary(model1412))$coefficients[,2]-(summary(model1412b))$coefficients[,2],3)
round((summary(model1413))$coefficients[,2]-(summary(model1413b))$coefficients[,2],3)
round((summary(model1414))$coefficients[,2]-(summary(model1414b))$coefficients[,2],3)
round((summary(model1415))$coefficients[,2]-(summary(model1415b))$coefficients[,2],3)


stargazer034<-htmlreg(list(model1409,model1410,model1411,model1412,model1413,model1414,model1415))

write(stargazer034,"results/Table A34.html")

#####

# TABLE A35
# split the data by country
# country = 1 is Denmark; country = 0 is the Netherlands

NLDK<-as.data.frame(as.matrix(cbind(RightDif,SelfLR,FEMAL,AGEL,BAMAL,Seat2,Spc,LvlReg,LvlNat,Country,ID1,ID2,BB,SB,BC,CSp,CS1,CS2,BI,NI,BN)))

modelnl21<-lmer(RightDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl22<-lmer(RightDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl23<-lmer(RightDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl24<-lmer(RightDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl25<-lmer(RightDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl26<-lmer(RightDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl27<-lmer(RightDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl28<-lmer(RightDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.

modelnl21b<-lmer(RightDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl22b<-lmer(RightDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl23b<-lmer(RightDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl24b<-lmer(RightDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl25b<-lmer(RightDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl26b<-lmer(RightDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl27b<-lmer(RightDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl28b<-lmer(RightDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])

round((summary(modelnl21))$coefficients[,1]-(summary(modelnl21b))$coefficients[,1],3)
round((summary(modelnl22))$coefficients[,1]-(summary(modelnl22b))$coefficients[,1],3)
round((summary(modelnl23))$coefficients[,1]-(summary(modelnl23b))$coefficients[,1],3)
round((summary(modelnl24))$coefficients[,1]-(summary(modelnl24b))$coefficients[,1],3)
round((summary(modelnl25))$coefficients[,1]-(summary(modelnl25b))$coefficients[,1],3)
round((summary(modelnl26))$coefficients[,1]-(summary(modelnl26b))$coefficients[,1],3)
round((summary(modelnl27))$coefficients[,1]-(summary(modelnl27b))$coefficients[,1],3)
round((summary(modelnl28))$coefficients[,1]-(summary(modelnl28b))$coefficients[,1],3)

round((summary(modelnl21))$coefficients[,2]-(summary(modelnl21b))$coefficients[,2],3)
round((summary(modelnl22))$coefficients[,2]-(summary(modelnl22b))$coefficients[,2],3)
round((summary(modelnl23))$coefficients[,2]-(summary(modelnl23b))$coefficients[,2],3)
round((summary(modelnl24))$coefficients[,2]-(summary(modelnl24b))$coefficients[,2],3)
round((summary(modelnl25))$coefficients[,2]-(summary(modelnl25b))$coefficients[,2],3)
round((summary(modelnl26))$coefficients[,2]-(summary(modelnl26b))$coefficients[,2],3)
round((summary(modelnl27))$coefficients[,2]-(summary(modelnl27b))$coefficients[,2],3)
round((summary(modelnl28))$coefficients[,2]-(summary(modelnl28b))$coefficients[,2],3)

stargazernl1<-htmlreg(list(modelnl21,modelnl22,modelnl23,modelnl24,modelnl25,modelnl26,modelnl27,modelnl28))

write(stargazernl1,"results/Table A35.html")

# TABLE A36

modelnl29<-lmer(RightDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl30<-lmer(RightDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl31<-lmer(RightDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl32<-lmer(RightDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl33<-lmer(RightDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl34<-lmer(RightDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl35<-lmer(RightDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


modelnl29b<-lmer(RightDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl30b<-lmer(RightDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl31b<-lmer(RightDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl32b<-lmer(RightDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl33b<-lmer(RightDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl34b<-lmer(RightDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl35b<-lmer(RightDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==0,])

round((summary(modelnl29))$coefficients[,1]-(summary(modelnl29b))$coefficients[,1],3)
round((summary(modelnl30))$coefficients[,1]-(summary(modelnl30b))$coefficients[,1],3)
round((summary(modelnl31))$coefficients[,1]-(summary(modelnl31b))$coefficients[,1],3)
round((summary(modelnl32))$coefficients[,1]-(summary(modelnl32b))$coefficients[,1],3)
round((summary(modelnl33))$coefficients[,1]-(summary(modelnl33b))$coefficients[,1],3)
round((summary(modelnl34))$coefficients[,1]-(summary(modelnl34b))$coefficients[,1],3)
round((summary(modelnl35))$coefficients[,1]-(summary(modelnl35b))$coefficients[,1],3)

round((summary(modelnl29))$coefficients[,2]-(summary(modelnl29b))$coefficients[,2],3)
round((summary(modelnl30))$coefficients[,2]-(summary(modelnl30b))$coefficients[,2],3)
round((summary(modelnl31))$coefficients[,2]-(summary(modelnl31b))$coefficients[,2],3)
round((summary(modelnl32))$coefficients[,2]-(summary(modelnl32b))$coefficients[,2],3)
round((summary(modelnl33))$coefficients[,2]-(summary(modelnl33b))$coefficients[,2],3)
round((summary(modelnl34))$coefficients[,2]-(summary(modelnl34b))$coefficients[,2],3)
round((summary(modelnl35))$coefficients[,2]-(summary(modelnl35b))$coefficients[,2],3)

stargazernl2<-htmlreg(list(modelnl29,modelnl30,modelnl31,modelnl32,modelnl33,modelnl34,modelnl35))

write(stargazernl2,"results/Table A36.html")

#####

# TABLE A37

modeldk21<-lmer(RightDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk22<-lmer(RightDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk23<-lmer(RightDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk24<-lmer(RightDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk25<-lmer(RightDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk26<-lmer(RightDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk27<-lmer(RightDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk28<-lmer(RightDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


modeldk21b<-lmer(RightDif~SelfLR+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk22b<-lmer(RightDif~SelfLR+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk23b<-lmer(RightDif~SelfLR+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk24b<-lmer(RightDif~SelfLR+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk25b<-lmer(RightDif~SelfLR+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk26b<-lmer(RightDif~SelfLR+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk27b<-lmer(RightDif~SelfLR+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk28b<-lmer(RightDif~SelfLR+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])

round((summary(modeldk21))$coefficients[,1]-(summary(modeldk21b))$coefficients[,1],3)
round((summary(modeldk22))$coefficients[,1]-(summary(modeldk22b))$coefficients[,1],3)
round((summary(modeldk23))$coefficients[,1]-(summary(modeldk23b))$coefficients[,1],3)
round((summary(modeldk24))$coefficients[,1]-(summary(modeldk24b))$coefficients[,1],3)
round((summary(modeldk25))$coefficients[,1]-(summary(modeldk25b))$coefficients[,1],3)
round((summary(modeldk26))$coefficients[,1]-(summary(modeldk26b))$coefficients[,1],3)
round((summary(modeldk27))$coefficients[,1]-(summary(modeldk27b))$coefficients[,1],3)
round((summary(modeldk28))$coefficients[,1]-(summary(modeldk28b))$coefficients[,1],3)

round((summary(modeldk21))$coefficients[,2]-(summary(modeldk21b))$coefficients[,2],3)
round((summary(modeldk22))$coefficients[,2]-(summary(modeldk22b))$coefficients[,2],3)
round((summary(modeldk23))$coefficients[,2]-(summary(modeldk23b))$coefficients[,2],3)
round((summary(modeldk24))$coefficients[,2]-(summary(modeldk24b))$coefficients[,2],3)
round((summary(modeldk25))$coefficients[,2]-(summary(modeldk25b))$coefficients[,2],3)
round((summary(modeldk26))$coefficients[,2]-(summary(modeldk26b))$coefficients[,2],3)
round((summary(modeldk27))$coefficients[,2]-(summary(modeldk27b))$coefficients[,2],3)
round((summary(modeldk28))$coefficients[,2]-(summary(modeldk28b))$coefficients[,2],3)

stargazerdk1<-htmlreg(list(modeldk21,modeldk22,modeldk23,modeldk24,modeldk25,modeldk26,modeldk27,modeldk28))

write(stargazerdk1,"results/Table A37.html")

# TABLE A38

modeldk29<-lmer(RightDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk30<-lmer(RightDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk31<-lmer(RightDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk32<-lmer(RightDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk33<-lmer(RightDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk34<-lmer(RightDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk35<-lmer(RightDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,])

# These models report "boundary (singular) fit: see help('isSingular')". As above we ran the model without respondent random effects and found that the results are identical to three decimal places.


modeldk29b<-lmer(RightDif~SelfLR+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk30b<-lmer(RightDif~SelfLR+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk31b<-lmer(RightDif~SelfLR+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk32b<-lmer(RightDif~SelfLR+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk33b<-lmer(RightDif~SelfLR+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk34b<-lmer(RightDif~SelfLR+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])
modeldk35b<-lmer(RightDif~SelfLR+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID2),data=NLDK[NLDK$Country==1,])

round((summary(modeldk29))$coefficients[,1]-(summary(modeldk29b))$coefficients[,1],3)
round((summary(modeldk30))$coefficients[,1]-(summary(modeldk30b))$coefficients[,1],3)
round((summary(modeldk31))$coefficients[,1]-(summary(modeldk31b))$coefficients[,1],3)
round((summary(modeldk32))$coefficients[,1]-(summary(modeldk32b))$coefficients[,1],3)
round((summary(modeldk33))$coefficients[,1]-(summary(modeldk33b))$coefficients[,1],3)
round((summary(modeldk34))$coefficients[,1]-(summary(modeldk34b))$coefficients[,1],3)
round((summary(modeldk35))$coefficients[,1]-(summary(modeldk35b))$coefficients[,1],3)

round((summary(modeldk29))$coefficients[,2]-(summary(modeldk29b))$coefficients[,2],3)
round((summary(modeldk30))$coefficients[,2]-(summary(modeldk30b))$coefficients[,2],3)
round((summary(modeldk31))$coefficients[,2]-(summary(modeldk31b))$coefficients[,2],3)
round((summary(modeldk32))$coefficients[,2]-(summary(modeldk32b))$coefficients[,2],3)
round((summary(modeldk33))$coefficients[,2]-(summary(modeldk33b))$coefficients[,2],3)
round((summary(modeldk34))$coefficients[,2]-(summary(modeldk34b))$coefficients[,2],3)
round((summary(modeldk35))$coefficients[,2]-(summary(modeldk35b))$coefficients[,2],3)

stargazerdk2<-htmlreg(list(modeldk29,modeldk30,modeldk31,modeldk32,modeldk33,modeldk34,modeldk35))

write(stargazerdk2,"results/Table A38.html")

# TABLE A39

model0401<-glmer(family=binomial(link = "logit"),Dic~Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0402<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0403<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0404<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0405<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0406<-glmer(family=binomial(link = "logit"),Dic~Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0407<-glmer(family=binomial(link = "logit"),Dic~Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0408<-glmer(family=binomial(link = "logit"),Dic~Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer011<-htmlreg(list(model0401,model0402,model0403,model0404,model0405,model0406,model0407,model0408))

write(stargazer011,"results/Table A39.html")

# TABLE A40

model0409<-glmer(family=binomial(link = "logit"),Dic~Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0410<-glmer(family=binomial(link = "logit"),Dic~Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0411<-glmer(family=binomial(link = "logit"),Dic~Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0412<-glmer(family=binomial(link = "logit"),Dic~Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0413<-glmer(family=binomial(link = "logit"),Dic~Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0414<-glmer(family=binomial(link = "logit"),Dic~Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0415<-glmer(family=binomial(link = "logit"),Dic~Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer012<-htmlreg(list(model0409,model0410,model0411,model0412,model0413,model0414,model0415))

write(stargazer012,"results/Table A40.html")

# TABLE A41

model0501<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0502<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0503<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0504<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0505<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0506<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0507<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0508<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer013<-htmlreg(list(model0501,model0502,model0503,model0504,model0505,model0506,model0507,model0508))

write(stargazer013,"results/Table A41.html")

# TABLE A42
model0509<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0510<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0511<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0512<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0513<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0514<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0515<-glmer(family=binomial(link = "logit"),Dic~Ext+Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer014<-htmlreg(list(model0509,model0510,model0511,model0512,model0513,model0514,model0515))

write(stargazer014,"results/Table A42.html")

# TABLE A43

model0601<-glmer(family=binomial(link = "logit"),Dic~FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0602<-glmer(family=binomial(link = "logit"),Dic~BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0603<-glmer(family=binomial(link = "logit"),Dic~BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0604<-glmer(family=binomial(link = "logit"),Dic~BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0605<-glmer(family=binomial(link = "logit"),Dic~BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0606<-glmer(family=binomial(link = "logit"),Dic~CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0607<-glmer(family=binomial(link = "logit"),Dic~CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0608<-glmer(family=binomial(link = "logit"),Dic~CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer015<-htmlreg(list(model0601,model0602,model0603,model0604,model0605,model0606,model0607,model0608))

write(stargazer015,"results/Table A43.html")

# TABLE A44

model0609<-glmer(family=binomial(link = "logit"),Dic~BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0610<-glmer(family=binomial(link = "logit"),Dic~SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0611<-glmer(family=binomial(link = "logit"),Dic~BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0612<-glmer(family=binomial(link = "logit"),Dic~BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0613<-glmer(family=binomial(link = "logit"),Dic~NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0614<-glmer(family=binomial(link = "logit"),Dic~BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model0615<-glmer(family=binomial(link = "logit"),Dic~BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer016<-htmlreg(list(model0609,model0610,model0611,model0612,model0613,model0614,model0615))

write(stargazer016,"results/Table A44.html")

# TABLE A45

model1401<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1402<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1403<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1404<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1405<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1406<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1407<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1408<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer015<-htmlreg(list(model1401,model1402,model1403,model1404,model1405,model1406,model1407,model1408))

write(stargazer015,"results/Table A45.html")

# TABLE A46

model1409<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1410<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1411<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1412<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1413<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1414<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1415<-glmer(family=binomial(link = "logit"),Dwc~Cwn1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer016<-htmlreg(list(model1409,model1410,model1411,model1412,model1413,model1414,model1415))

write(stargazer016,"results/Table A46.html")

# TABLE A47

model1501<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1502<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1503<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1504<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1505<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1506<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1507<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1508<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer017<-htmlreg(list(model1501,model1502,model1503,model1504,model1505,model1506,model1507,model1508))

write(stargazer017,"results/Table A47.html")

# TABLE A48

model1509<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1510<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1511<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1512<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1513<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1514<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1515<-glmer(family=binomial(link = "logit"),Dic~Pri+Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer018<-htmlreg(list(model1509,model1510,model1511,model1512,model1513,model1514,model1515))

write(stargazer018,"results/Table A48.html")

# TABLE A49

model1601<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1602<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1603<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1604<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1605<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1606<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1607<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1608<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer019<-htmlreg(list(model1601,model1602,model1603,model1604,model1605,model1606,model1607,model1608))

write(stargazer019,"results/Table A49.html")

# TABLE A50

model1609<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1610<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1611<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1612<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1613<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1614<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))
model1615<-glmer(family=binomial(link = "logit"),Dic~EXP5+Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+Country+(1|ID1)+(1|ID2))

stargazer020<-htmlreg(list(model1609,model1610,model1611,model1612,model1613,model1614,model1615))

write(stargazer020,"results/Table A50.html")

#######
# TABLE A51
# split the data by country
# country = 1 is Denmark; country = 0 is the Netherlands

NLDK<-as.data.frame(as.matrix(cbind(Dic,Con,Con1,FEMAL,AGEL,BAMAL,Seat2,Spc,LvlReg,LvlNat,Country,ID1,ID2,BB,SB,BC,CSp,CS1,CS2,BI,NI,BN)))

modelnl41<-glmer(family=binomial(link = "logit"),Dic~Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl42<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl43<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl44<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl45<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl46<-glmer(family=binomial(link = "logit"),Dic~Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl47<-glmer(family=binomial(link = "logit"),Dic~Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl48<-glmer(family=binomial(link = "logit"),Dic~Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])

stargazernl01<-htmlreg(list(modelnl41,modelnl42,modelnl43,modelnl44,modelnl45,modelnl46,modelnl47,modelnl48))

write(stargazernl01,"results/Table A51.html")

# TABLE A52

modelnl49<-glmer(family=binomial(link = "logit"),Dic~Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl50<-glmer(family=binomial(link = "logit"),Dic~Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl51<-glmer(family=binomial(link = "logit"),Dic~Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl52<-glmer(family=binomial(link = "logit"),Dic~Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl53<-glmer(family=binomial(link = "logit"),Dic~Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl54<-glmer(family=binomial(link = "logit"),Dic~Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])
modelnl55<-glmer(family=binomial(link = "logit"),Dic~Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==0,])

stargazernl02<-htmlreg(list(modelnl49,modelnl50,modelnl51,modelnl52,modelnl53,modelnl54,modelnl55))

write(stargazernl02,"results/Table A52.html")

# TABLE A53
# I set the number of points per axis for evaluating the adaptive Gauss-Hermite approximation to the log-likelihood (nAGQ) to zero to deal with convergence issues

modelDK41<-glmer(family=binomial(link = "logit"),Dic~Con1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK42<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0) ##
modelDK43<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0) ##
modelDK44<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK45<-glmer(family=binomial(link = "logit"),Dic~Con1+BC+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0) ##
modelDK46<-glmer(family=binomial(link = "logit"),Dic~Con1+CSp+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK47<-glmer(family=binomial(link = "logit"),Dic~Con1+CS1+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK48<-glmer(family=binomial(link = "logit"),Dic~Con1+CS2+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)

stargazerDK01<-htmlreg(list(modelDK41,modelDK42,modelDK43,modelDK44,modelDK45,modelDK46,modelDK47,modelDK48))

write(stargazerDK01,"results/Table A53.html")

# TABLE A54
# nAGQ=0 to deal with convergence issues
modelDK49<-glmer(family=binomial(link = "logit"),Dic~Con1+BB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK50<-glmer(family=binomial(link = "logit"),Dic~Con1+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK51<-glmer(family=binomial(link = "logit"),Dic~Con1+BB+SB+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK52<-glmer(family=binomial(link = "logit"),Dic~Con1+BI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK53<-glmer(family=binomial(link = "logit"),Dic~Con1+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)
modelDK54<-glmer(family=binomial(link = "logit"),Dic~Con1+BI+NI+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0) ##
modelDK55<-glmer(family=binomial(link = "logit"),Dic~Con1+BN+FEMAL+AGEL+BAMAL+Seat2+Spc+LvlReg+LvlNat+(1|ID1)+(1|ID2),data=NLDK[NLDK$Country==1,],nAGQ=0)

stargazerDK02<-htmlreg(list(modelDK49,modelDK50,modelDK51,modelDK52,modelDK53,modelDK54,modelDK55))

write(stargazerDK02,"results/Table A54.html")

# TABLE A55
# this is split per issue, so we had to unstack a number of variables

model1401<-(lm(RightDif[1:2994]~SelfLR[1:2994]+PriASYL+FEMA+Age[1:2994]+BAMA+Seats[1:2994]+SpcASYL+LvlREG+LvlNAT+Countr))
model1402<-(lm(RightDif[2994+c(1:2994)]~SelfLR[2994+c(1:2994)]+PriCCTV+FEMA+Age[1:2994]+BAMA+Seats[1:2994]+SpcCCTV+LvlREG+LvlNAT+Countr))
model1403<-(lm(RightDif[2994*2+c(1:2994)]~SelfLR[2994*2+c(1:2994)]+PriCULT+FEMA+Age[1:2994]+BAMA+Seats[1:2994]+SpcCULT+LvlREG+LvlNAT+Countr))
model1404<-(lm(RightDif[2994*3+c(1:2994)]~SelfLR[2994*3+c(1:2994)]+PriCARE+FEMA+Age[1:2994]+BAMA+Seats[1:2994]+SpcCARE+LvlREG+LvlNAT+Countr))
model1405<-(lm(RightDif[2994*4+c(1:2994)]~SelfLR[2994*4+c(1:2994)]+PriROAD+FEMA+Age[1:2994]+BAMA+Seats[1:2994]+SpcROAD+LvlREG+LvlNAT+Countr))

stargazer037<-htmlreg(list(model1401,model1402,model1403,model1404,model1405))

write(stargazer037,"results/Table A55.html")

########################################################################################
# PART 5: MEDIATION ANALYSES
# NOTE: THE MEDIATIONS ARE SIMULATED AND CAN DIFFER PER ANALYSIS.

# TABLE A56: Model 373
# generate a data set with the relevant variables
datamed01<-as.data.frame(cbind(Abs,Ext,Con,RightDif,SelfLR,BC,CSp,Spc,FEMAL,AGEL,BAMAL,Seat2,LvlReg,LvlNat,Country,ID1,ID2))
datamed02<-datamed01[rowSums(is.na(datamed01))==0,]
datamed03<-as.list(datamed02)

model01<-(lmer(RightDif~SelfLR+BC+CSp+Spc+FEMAL+AGEL+BAMAL+Seat2+LvlReg+LvlNat+Country+(1|ID2),data=datamed03))

stargazerMD01<-htmlreg(list(model01))

write(stargazerMD01,"results/Table A56 - 373 - 1.html")

medfit01<-(lmer(SelfLR~BC+(1|ID2),data=datamed03))

stargazerMD02<-htmlreg(list(medfit01))

write(stargazerMD02,"results/Table A56 - 373 - 2.html")

# This is the mediation, that is copied manually into the table 
med.out01 <- mediate(medfit01, model01, treat = "BC", mediator = "SelfLR",
                  sims = 10000)                 
summary(med.out01)

# TABLE A56: Model 374

medfit02<-(lmer(SelfLR~CSp+(1|ID2),data=datamed03))

stargazerMD03<-htmlreg(list(medfit02))

write(stargazerMD03,"results/Table A56 - 374 - 2.html")

# This is the mediation, that is copied manually into the table
med.out02 <- mediate(medfit02, model01, treat = "CSp", mediator = "SelfLR",
                  sims = 10000)                
summary(med.out02)

# TABLE A56: Model 375

medfit03<-(lmer(BC~SelfLR+(1|ID2),data=datamed03))

stargazerMD04<-htmlreg(list(medfit03))

write(stargazerMD04,"results/Table A56 - 375 - 2.html")

# This is the mediation, that is copied manually into the table 
med.out03 <- mediate(medfit03, model01, treat = "BC", mediator = "SelfLR",
                  sims = 10000)                 
summary(med.out03)

# TABLE A56: Model 376
medfit04<-(lmer(CSp~SelfLR+(1|ID2),data=datamed03))

stargazerMD05<-htmlreg(list(medfit04))

write(stargazerMD05,"results/Table A56 - 376 - 2.html")

# This is the mediation, that is copied manually into the table 
med.out04 <- mediate(medfit04, model01, treat = "CSp", mediator = "SelfLR",
                  sims = 10000)                
summary(med.out04)
 
# TABLE A57: Model A377

model02<-(lmer(Abs~Con+BC+CSp+Spc+FEMAL+AGEL+BAMAL+Seat2+LvlReg+LvlNat+Country+(1|ID2),data=datamed03))

stargazerMD06<-htmlreg(list(model02))

write(stargazerMD06,"results/Table A57 - 377 - 1.html")

medfit05<-(lmer(Con~BC+(1|ID2),data=datamed03))

stargazerMD07<-htmlreg(list(medfit05))

write(stargazerMD07,"results/Table A57 - 377 - 2.html")
       
# This is the mediation, that is copied manually into the table    
med.out05 <- mediate(medfit05, model02, treat = "BC", mediator = "Con",
                  sims = 10000)                 
summary(med.out05)

# TABLE A57: Model A378

medfit06<-(lmer(Con~CSp+(1|ID2),data=datamed03))

stargazerMD08<-htmlreg(list(medfit06))

write(stargazerMD08,"results/Table A57 - 378 - 2.html")
    
# This is the mediation, that is copied manually into the table
med.out06 <- mediate(medfit06, model02, treat = "CSp", mediator = "Con",
                  sims = 10000)  
summary(med.out06)     

# TABLE A57: Model 379

medfit07<-(lmer(BC~Con+(1|ID2),data=datamed03))

stargazerMD09<-htmlreg(list(medfit07))

write(stargazerMD09,"results/Table A57 - 379 - 2.html")

# This is the mediation, that is copied manually into the table
med.out07 <- mediate(medfit07, model02, treat = "BC", mediator = "Con",
                  sims = 10000)                 
summary(med.out07) 

# TABLE A57: Model 380    
medfit08<-(lmer(CSp~Con+(1|ID2),data=datamed03))

stargazerMD10<-htmlreg(list(medfit08))

write(stargazerMD10,"results/Table A57 - 380 - 2.html")

# This is the mediation, that is copied manually into the table
med.out08 <- mediate(medfit08, model02, treat = "CSp", mediator = "Con",
                  sims = 10000)  
summary(med.out08)      

###########
# THE END #
########### 